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ABSTRACT 


A  numerical  method  is  presented  foi  the  analysis  of  steady-state 
wave  propagation  problems  in  linearly  elastic  or  viscoelastic  media  of 
infinite  extent.  Plane  and  axisymmetric  geometries  are  considered  which 
consist  of  a  finite  irregular  region  joined  to  semi-infinite  layered 
regions.  By  this  method,  torsional  and  vertical  vibrations  of  circular 
footings  on,  or  embedded  in,  homogeneous  and  inhomogeneous  soil  layers 
over  rock  are  studied. 

The  irregular  region  is  discretized  by  compatible  finite  elements, 
while  the  semi-infinite  layered  regions  are  discretized  by  subdividing 
the  layers  into  thin  sublayers  and  by  assuming  that  within  each  sublayer 
the  displacements  vary  linearly  in  the  direction  normal  to  the  layers. 

In  the  direction  parallel  to  the  layers,  the  displacements  are  expanded 
into  a  finite  number  of  plane  or  axisymmetric  propagating  and  decaying 
wave  modes  which  are  determined  by  the  solution  of  algebraic  eigenvalue 
problems.  Dynamic  stiffness  matrices  are  developed  which  uniquely  relate 
nodal  forces  to  simultaneous  noial  displacements  at  the  boundary  between 
the  irregular  and  the  layered  regions  and  thus  represent  the  dynamic 
response  of  the  semi-infinite  layered  regions.  The  dynamic  stiffness  of 
the  combined  regions  is  computed  by  the  direct  stiffness  method.  The 
equations  of  motion  are  derived  from  the  principle  of  virtual  work  which 
is  formulated  to  facilitate  the  use  of  complex  variables  in  the  analysis 
of  the  harmonic  motion. 

Solutions  obtained  by  this  numerical  method  show  good  agreement 
with  known  analytical  solutions.  The  response  of  circular  footings, 
which  are  supported  by  soil  layers  over  rock  and  are  excited  to  vibrate 
in  the  torsional  and  vertical  mode,  are  investigated;  the  effects  of  the 
thickness  of  the  supporting  layer,  embedment  of  the  footing  and  increas¬ 
ing  shear  modulus  with  depth  are  studied.  The  screening  effect  of 
trenches  on  horizontally  polarized  shear  waves  is  explored  by  varying 
the  trench  depth  and  the  frequency  of  the  wave  motion. 
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FOREWORD 


The  research  described  in  this  report  was  performed  -under  Contract 
?5c..  DACA39-70-C-0023,  "Dynamic  Analysis  of  Embedded  Footings  on  Layered 
Subsoil,”  between  the  U.  S.  Army  Engineer  Waterways  Experiment  Station 
(WES)  and  The  Regents  of  the  University  of  California  during  the  period 
1  -June  1970  to  30  -Tune  1972. 

The  objective  of  the  research,  which  was  begun  in  June  1970,  was 
to  compute  st«  ady  state  displacements  generated  within  a  layered  soil 
system  by  forced  vibrations  of  a  rigid  circular  footing  located  on,  or 
embedded  in,  this  sytem.  The  stated  objective  has  been  achieved  and  in 
addition  a  nev  very  general  method  has  been  developed  for  steady  state 
analysis  of  ir finite  plane  and  axisymmetric  structures.  The  major  part 
of  the  research  described  was  performed  by  Mr.  Gunter  Waas  as  part  of 
his  ioctoral  research  under  the  supervision  of  the  principal  investiga¬ 
tor,  Dr.  John  bysmer,  Associate  Professor  of  Civil  Engineering. 

The  oroject  was  administered  by  the  Office  of  Research  Services 
of  the  College  of  Engineering.  The  contracting  officer  for  WES  was 
Mr.  John  J.  Kirschenbaun,  Jr.,  represented  by  Dr.  Lyman  W.  Heller  of 
the  Soils  and  Pavements  Laboratory.  During  the  period  of  this  work, 

Mr.  J.  P.  Sale  was  Chief  of  the  Soils  and  Pavements  Laboratory,  Mr. 

R.  G.  Ahlvin  was  Assistant  Chief  of  the  Soils  and  Pavements  Laboratory, 
and  Mr.  R.  W.  Cunny  was  Chief  of  the  Soil  Dynamics  Branch. 

CCL  Ernest  Peixotto,  CE,  was  the  Director  of  the  WES  and 
Mr.  F.  R.  Brovn  was  the  Technical  Director. 
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1.  INTRODUCTION 


1.1  Foundation  Vibration  Problems 

The  designer  of  foundations  for  machinery  must  be  able  to  predict 
whether  or  not  the  response  of  a  foundation  to  dynamic  loads  from  the 
machine  will  meet  certain  design  criteria  and  if  the  vibrations  trans¬ 
mitted  through  the  ground  will  be  small  enough  as  not  to  affect  other 
machines,  vibration  sensitive  equipment  and  humans  in  the  neighborhood. 
The  prediction  of  the  vibration  amplitudes  involves  the  determination 
of  the  dynamic  loads  and  soil  properties,  the  selection  of  a  mathematical 
model  and  its  analysis. 

Discussions  and  data  on  dynamic  loads  from  machines  are  presented 
in  Refs.  [9,  27].  Periodic  loads  caused  by  combustion  engines,  com¬ 
pressors  and  turbins  are  to  be  distinguished  from  transient  loads  caused 
by  forge  hammers  and  testing  machines.  If  the  soil- foundation  system 
behaves  linearly  under  dynamic  loading,  the  response  to  periodic  loads 
can  be  obtained  by  superposition  of  the  harmonic  response  at  different 
frequencies.  Trans1' o-  t  loads  may  be  t  .•eate-’  similarly  by  employing  the 
Fourier  transfv rmation,  possibly  in  a  discretized  form. 

The  determination  of  the  dynamic  soil  properties  is  discussed 
briefly  in  Chapter  2.  The  dynamic  stress-strain  behavior  of  soils  is 
approximately  linear  for  small  strain  amplitudes  and  can  therefore  be 
represented  by  elastic  moduli.  To  incorporate  material  damping  in 
harmon"  '.on  complex  moduli  may  be  employed. 

iometry  involved  is  often  quite  complicated.  It  may  consist 
of  the  foundation,  several  soil  layers  with  different  material  proper¬ 
ties,  rock  at  some  depth  and  adjacent  foundations  and  structures.  To 
make  the  problem  amenable  to  analysis  the  geometry  must  therefore  be 
strongly  idealized. 

Many  problems  which  are  actually  three-dimensional  in  space  can 
be  reduced  with  some  justification  to  either  plane  or  axisymmetric 
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problems-  For  instance,  the  soil-foundation  interaction  of  a  turbine 
pedestal,  which  is  ten  times  as  long  as  it  is  wide,  may  be  treated  as 
a  plane  problem  if  the  motion  of  the  foundation  slab  in  a  plane  per¬ 
pendicular  to  the  turbine  axis  is  studied.  Thus  a  cross-section  of  the 
soil-foundation  system  is  considered  which  is  plane  and  of  unit  thick¬ 
ness.  On  the  othei  hand,  the  vertical  vibration  o.t  u  machine  foundation 
which  is  square  in  plan  may  be  treated  as  an  axisymmecric  problem.  For 
this  purpose  the  square  footing  is  approximated  by  a  circular  footing 
of  equal  area  in  plan. 

Many  problems  in  soil  dynamics  can  be  characterized  as  wave 
propagation  problems,  because  the  motion  in  the  soil  is  dominated  by 
propagating  waves.  The  waves  are  generated  at  some  source,  for  instance 
by  the  vibration  of  a  machine  footing,  and  from  there  they  propagate 
into  the  underlying  medium  which  is  unbounded  except  at  the  surface  of 
the  ground.  The  effect  of  the  propagating  waves  is  to  disperse  energy 
into  the  infinite  medium  and  at  the  same  time  to  carry  disturbances  a 
long  distance  away  from  t  le  source.  Due  to  energy  dissipation  by  the 
waves,  the  footing  response  to  dynamic  loads  is  damped  and  the  displace¬ 
ment  amplitudes  of  the  footing  are  considerably  reduced,  especially 
at  frequencies  corresponding  to  oximum  amplitudes. 

The  depth  profile  at  a  construction  site  commonly  consists  of 
Horizontally  layered  soil  deposits  underlain  by  rock.  The  lower  soi2 
layers  are  usually  stiffer  than  the  upper  layers,  because  the  overburden 
and  thus  the  confining  or  consolidation  pressure  increases  with  depth. 
Waves  generated  somewhere  close  to  the  surface  of  the  ground  are 
therefore  reflected  and  refracted  at  the  layer  interfaces.  Since  the 
soil  layers  generally  become  stiffer  with  depth,  the  pattern  of  reflec¬ 
tions  and  refractions  causes  the  waves  to  propagate  mainly  in  the 
horizontal  direction  and  most  of  the  dispersed  energy  is  transmitted 
in  the  softer  layers  close  to  the  surface.  Even  in  the  case  of  a 
vertically  vibrating  circular  footing  which  is  supported  by  a  homo¬ 
geneous  elastic  half  space.  Fig.  1,  67  percent  cf  the  energy  is  dis¬ 
persed  by  Rayleigh  waves  which  travel  along  the  surface  of  the  half 
space  f28].  Therefore  an  even  higher  percentage  of  the  energy  will  be 
dispersed  in  the  horizontal  direction  if  the  medium  is  stratified  as 


? 


described  above.  If  the  rock  underlying  the  soil  deposits  is  much 
harder  than  the  soil,  as  is  usually  the  case,  it  may  be  assumed  for 
practical  purposes  that  the  rock  is  completely  rigid  and  that  the  waves 
are  totally  reflected  at  the  surface  of  the  rock. 


1.2  Existing  Methods  of  Analysis 


The  dynamic  analysis  of  machine  foundations  is  often  based  on  the 

X 

mathematical  model  shown  in  Fig.  1.  A  circular  footing  is  supported  by 
a  homogeneous,  isotropic,  linearly  elastic,  semi- inf inite  "half  space. 
This  model  was  first  studied  by  Reissner  [39,  40,  41]  in  the  1930’s. 

He  derived  analytical  solutions  for  the  forced  vertical  and  torsional 
vibra. ion  of  the  footing.  Subsequently,  the  model  has  been  studied  by 
some  20  researchers  who  have  directed  almost  their  entire  attention  to 
the  dynamic  compliances  of  the  vertical,  torsional,  rocking  ahd  hori¬ 
zontal  vibration  of  the  circular  footing.  A  recently  published  paper 
by  Luco  and  Westmann  [23]  presents  a  summary  of  the  various  .solutions 
and  numerical  results  for  the  compliances  over  a  large  frequency  range. 
The  paper  also  contains  some  information  on  the  stress  distribution 
under  the  footing  and  on  the  displacements  in  the  far  field. 

A  mere  satisfactory  model  for  many  foundation  vibration  problems 
where  rock  is  situated  at  shallow  depth  is  shown  in  Fig.  2.  Here,  a 
circular  footing  is  supported  by  a  homogeneous  elastic  layer  over  a 
rigid  base.  This  model  is  more  difficult  to  deal  with  by  analytical 
methods  than  the  half  space  model  because  additional  boundary  conditions 
have  to  be  satisfied  at  the  rigid  base.  The  forced  torsional  motion 
of  the  footing  was  studied  by  Arnold  et  al.  [6],  Bycroft  [10]  and 
Awojobi  [8].  The  forced  vertical  motion  of  the  footing  was  investigated 
by  Bycroft  [10]  and  Warburtcn  [46].  No  methods  of  analysis  for  the 
rocking  and  sliding  motion  are  available. 

Froblems  which  include  layering  of  the  soil  or  embedment  of  the 
foundation,  two  significant  features  of  most  foundation  vibration 
problems,  have  not  been  solved  by  analytical  methods,  since  the  in¬ 
herent  boundary  conditions  are  too  complex. 


2 


I 

I 


I 


Numerical  methods  such  as  the  finite  difference  and  the  finite 
element  method  are  not  directly  applicable  to  steady-state  wave 
propagation  problems  in  an  infinite  medium,  because  only  a  finite 
number  of  nodal  points  can  be  considered  and  therefore  the  discrete 
model  is  confined  to  a  finite  region.  This  poses  the  problem  as  to 
how  wave  reflections  can  be  prevented  at  an  artificial  boundary,  whicii 
is  introduced  to  confine  the  model. 

Lysmer  and  Kuhlemeyer  [25]  developed  a  viscous  boundary  which 
absorbes  much  of  the  impinging  wave  energy.  At  the  boundary  they 
applied  viscous  forces  which  they  computed  from  one-dimensional  wave 
propagation  theory.  The  dashpots  can  be  designed  to  perfectly  absorb 
the  energy  contained  in  P-waves  (compressions.!  waves)  and  S-waves 
(shear  waves)  If  the  incident  angles  of  the  waves  are  known.  However, 
the  incident  angles  are  usually  unknown  and  therefore  must  be 
anticipated.  If  the  actual  incident  angles  differ  from  those  anticipated 
by  not  more  than  about  60  degrees,  the  dashpots  are  very  effective  [25]. 

Lysmer  and  Kuhlemeyer  obtained  good  results  with  discrete  models 
of  the  type  shown  in  Fig.  3  when  applied  to  the  analysis  of  steady-state 
vertical  vibrations  of  circular  footings  embedded  in  a  homogeneous 
elastic  half  space.  They  designed  the  dashpots  at  the  lower  horizontal 
boundary  for  perpendicular  incident  angles  and  the  dashpots  at  the 
vertical  boundary  "or  the  frequencv  dependent  incident  angles  of  the 
P  and  S-wave  components  that  form  a  Rayleigh  wave. 

Kuhlemeyer  [22]  extended  the  method  to  the  case  of  a  layered 
half  space.  He  obtained  good  results  in  the  low  frequency  range  fer 
case  .  in  which  the  elastic  moduli  of  the  layers  decreased  wi *'h  depth. 
However,  he  observed  that  in  general  it  is  difficult  to  predict  the 
incident  angles  of  P  and  S -waves  at  the  artificial  boundary,  because 
the  ?  and  S-waves  undergo  multiple  reflections  and  refractions  at 
the  layer  interfaces  and  at  the  free  surface.  To  complicate  the  matter, 
the  incident  angles  often  change  drastically  with  frequency. 

Ang  and  Newnark  [4]  developed  a  "transmitting  boundary"  at  about 
the  same  timp  that  Lysmer  and  Kuhlemeyer  developed  their  viscous 
boundary.  In  ssence,  both  boundaries  are  the  same.  They  are  based  on 
one-dimens ionai  wave  propagation  theory  and  require  some  advance 


knowledge  of  the  incident  angles.  Ang  and  Newmark  used  the  transmitting 
boundary  in  the  analysis  of  ground  shock  waves  caused  by  nuclear  blasts. 
They  demonstrated  the  effectiveness  of  the  boundary  in  cases  in  which 
the  wave  motion  was  mainly  one-dimensional,  but  they  presented  little 
evidence  of  the  boundary's  performance  in  two-dimensional  wave  propaga¬ 
tion  problems. 

However,  studies  by  Hadala  [IS]  show  that  the  transmitting  boundary 
is  useful  in  certain  two-dimensional  shock  wave  problems  in  layered 
semi-infinite  media,  when  only  a  very  short  time  period  after  the  blast 
is  considered.  The  degree  of  perfection  required  of  a  transmitting 
boundary  in  these  problems  is  less  string!...:  than  in  steady-state  wave 
propagation  problems. 

Another  method  of  analyzing  steady-state  wave  propagation  problems 
in  infinite  media  is  ,_o  employ  a  very  large  discrete  model  with  a  con¬ 
siderable  amount  of  damping.  Thus,  the  waves  generated  at  the  source 
are  strongly  attenuated  before  they  reach  the  artificial  boundary,  and 
after  reflection,  are  further  attenuated  so  that  they  are  almost 
dissipated  before  they  return  to  the  zone  of  interest.  This  approach 
is  physically  justified  if  the  chosen  amount  of  damping  is  reasonable 
and  ii  Lha  model  is  large  enough .  B-.t  ii.  pz'ctiee,  it  is  p-wbibitive’y 
expensive  in  computational  time  and  storage  requirements  despite  the 
dramatic  advances  in  the  development  of  digital  computers. 

A  new  discrete  model  which  is  particularly  suited  for  the  analysis 
of  foundation  vibration  problems  is  described  below. 

1  *  3  A  Finite  Dynamic  Model  for  Layered  Semi-Infinite  Media 

The  discrete  method  developed  in  this  dissertation  is  designed  to 
solve  plane  and  axisynanetric  wave  propagation  problems  in  semi-inf inite, 
linearly  viscoelastic  media  of  the  type  illustrated  in  Figs.  4  and  5. 

Fig.  4  shows  the  cross-section  of  a  typical  plane  problem  in 
which  the  geometry,  the  material  properties  and  the  external  loads  do 
not  change  in  the  direction  perpendicular  to  the  cross-section.  The 
region  1  is  joined  along  vertical  planes  to  the  regions  L  and  R  which 
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extend  infinitely  to  the  left  and  to  the  right,  respectively. 

The  regions  L  and  8  consist  of  horizontal  layers  which  are  welded 
together  at  their  interfaces  and  may  differ  in  their  material  properties 
and  thickness.  They  are  called  the  layered  or  semi-infinite  regions. 

The  region  I  contains  all  geometrical  irregularities  and  is  therefore 
called  the  irregular  region.  The  regions  L,  R  and  1  are  welded  to  a 
rigid  base  which  is  fixed  in  space. 

The  materials  involved  may  be  linearly  elastic  or  linearly  visco¬ 
elastic.  For  simplicity,  they  are  assumed  to  be  isotropic. 

All  external  forces  act  within  the  irregular  regions  and  vary 
harmonically  with  time.  They  generate  harmonic  motion  which  consists  of 
propagating  and  standing  waves.  The  propagating  waves  transmit  energy. 
This  energy  is  either  dissipated  due  to  material  damping  within  a 
finite  distance  from  its  source  or  is  propagated  outward  to  infinity 
if  the  layered  regions  are  ideally  elastic. 

Fig.  5  shows  a  typical  cross-section  of  the  analogous  axisymmetric 
problem  in  which  the  geometry,  the  material  properties  and  the  loads  do 
not  vary  in  the  angular  direction. 

Since  there  is  no  hope  tu  find  closed  form  solutions  to  problems 
as  complicated  as  this,  nuiueiical  methods  which  reduce,  the  difficulty 
by  discretization  must  be  resorted  to. 

The  finite  element  method  is  the  most  flexible  numerical  method 
available  for  solving  complex  boundary  value  problems  in  continuum 
mechanics,  because  it  is  easily  adaptable  to  complicated  geometries 
and  local  changes  in  material  properties.  Therefore  it  is  employed  to 
analyze  the  motion  in  the  irregular  region.  It  is  applied  to  the 
present  problem  in  the  displacement  formulation.  The  application  is 
outlined  in  Chapter  3. 

The  layered  regions,  however,  cannot  be  treated  by  the  finite 
element  method,  because  these  regions  are  of  infinite  extent  and  their 
discretization  by  finite  elements  would  result  in  an  infinite  number 
of  elements  and  degrees  of  freedom. 

Since  the  geometry  of  the  layered  region,  L  or  R,  does  not  vary 
in  the  horizontal  direction,  a  numerical  method  is  suggested  which 
reduces  the  analysis  of  harmonic  motion  in  the  layered  region  to  a 


discrete  problem  with  typically  20  to  60  degrees  of  freedom.  To  this 
end,  the  region  is  discretized  in  the  following  way.  The  natural 
layers  are  subdivided  into  thin  sublayers  as  shown  in  Fig.  6  and  it  is 
assumed  that  the  displacements  within  each  sublayer  vary  linearly  in 
the  vertical  direction.  But  in  the  horizontal  direction  the  displace¬ 
ments  are  required  to  satisfy  the  pertinent,  ordinary  differential 
equations  which  are  obtained  from  the  governing,  partial  differential 
equations  by  separating  the  variables.  This  method  is  similar  to 
Lysmer's  original  lumped  mass  method  for  the  analysis  of  Rayleigh 
waves  [26]. 

At  a  given  frequency,  the  free  motion  in  the  thus  discretized 
layered  region  consists  of  a  finite  number  of  wave  modes  which  are 
obtained  by  the  solution  of  an  algebraic  eigenvalue  problem.  These 
wave  modes  serve  as  shape  functions  for  expanding  the  displacements  in 
the  region  in  terms  of  mode  participation  factors. 

If  the  layered  region  is  separated  from  the  irregular  region, 
nodal  forces  have  to  be  applied  at  the  boundary  as  shown  in  Figs.  11  and 
12  in  order  to  preserve  dynamic  equilibrium.  These  nodal  forces  are 
derived  from  the  displacement  expansion  by  observing  the  strain- 
Jicpla..eme.il  and  s i r a x. . i  relalaons.  ine  noda.  forces  are 

uniquely  related  to  the  simultaneous  nodal  displacements  at  the  boundary 
through  a  dynamic  stiffness  matrix  which  represents  the  elastic  and 
the  viscous  response  of  the  semi-infir.ite  region.  The  dynamic  stiffness 
matrix  is  complex  symmetric  and  frequency  dependent. 

Once  thi  dynamic  stiffness  matrices  of  the  layered  regions  are 
established,  the  analysis  of  the  combined  regions  L,  I  and  R  follows 
the  usual  procedure  of  the  direct  stiffness  method.  The  analysis 
yields  the  displacements  and  stresses  in  the  irregular  and  the  layered 
regions. 

The  dynamic  stiffness  matrix  of  a  layered  region  is  developed  in 
Chapters  A  through  7  for  the  plane  and  the  axisymmetric  case.  In  each 
case  the  motions  in  the  plane  and  perpendicular  to  the  plane  of  the 
cross-section,  shown  in  Figs.  A  and  5,  are  treated  separately  since 
they  are  independent. 


In  the  layered  regions,  the  motions  in  the  plane  and  perpendicular 
to  the  plane  of  the  cross-section  consist  of  generalized  Rayleigh  and 
Love  waves,  respectively. 

A  generalized  Rayleigh  wave  in  the  plane  case  is  here  defined  by 
6  =  u(z)  exp  (iwt  -  ikx) 

X  (1.1) 

5  =  w(z)  exp  (iwt  -  ikx) 

z 

in  which  6  and  6  are  the  displacements  in  the  x  and  z-direction,  u(z) 
x  z  r 

and  w(z)  are  the  corresponding  vertical  mode  shapes,  w  is  the  frequency, 
k  is  called  the  wave  number  and  i  =  </-l.  The  wave  defined  by  Eq.  (1.1) 
is  called  a  generalized  Rayleigh  wave,  because  Rayleigh  waves  in 
layered  media  are  described  mathematically  in  the  same  way,  but  are 
usually  understood  as  having  real  wave  numbers  [16],  whereas  the  wave 
number  in  Eq.  (1.1)  may  be  real,  imaginary  or  complex. 

A  generalized  Love  wave  in  the  plane  case  is  defined  by 

<5y  =  v(z)  exp  (iwt  -  ikx)  (1.2) 

in  which  6^  is  the  displacement  in  the  y-direction  and  v(z)  is  the 
vertical  mode  shape.  The  wave  number  may  again  be  real,  imaginary  or 
crmplex,  while  Love  waves  in  layered  media  are  usually  understood  as 
having  real  wave  '-".rubers  [16]. 

Axisymmetric  generalized  Rayleigh  and  Love  waves  are  described 
similarly  in  Chapters  5  and  7.  Henceforth,  the  cases  in  which  the 
motion  occurs  either  in  the  plane  or  perpendicular  to  the  plane  of  the 
cross-section  are  called  the  Rayleigh  wave  case  or  the  Love  wave  case 
according  to  the  type  of  wave  generated. 
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2.  DYNAMIC  SOIL  PROPERTIES 


2.1  Introduction 

Successful  determination  of  the  response  of  soil-foundation 
systems  to  dynamic  loads  is  critically  dependent  on  the  incorporation 
of  representative  soil-properties  in  the  analysis.  Thus  considerable 
effort  has  been  directed  towards  the  determination  of  soil  properties 
in  recent  years.  Several  techniques  for  measuring  dynamic  soil  proper¬ 
ties  have  been  either  newly  developed  or  improved  and  a  considerable 
amount  of  data  has  been  collected  [42,  43]. 

The  stress-strain  relationships  of  most  soils  subjected  to 
symmetric  cyclic  loading  conditions  are  curvi-linear  as  shown  in  Fig.  7. 
The  secant  modulus,  which  is  determined  by  the  extreme  points  on  the 
hysteresis  loop,  and  the  damping  factor,  which  is  proportional  to  the 
area  inside  the  loop,  may  be  select'*-’  to  represent  the  hysteretic 
stress-strain  behavior  in  analyses.  Both  the  modulus  and  the  damping 
factor  depend  tc  some  degree  on  the  magnitude  of  the  strain  as  the 
axis  and  the  area  of  the  loop  vary  with  the  strain  amplitude. 


2.2  Experimental  Determination 


The  main  procedures  for  determining  the  moduli  and  damping 
characteristics  at  low  to  moderately  high  strain  levels  are  summarized 
below  [43]. 

a)  Direct  determination  of  stress-strain  relationships 

Hysteretic  stress-strain  curves  of  the  type  shown  in  Fig.  / 

can  be  obtained  in  the  laboratory  from  triaxial  compression 

tests,  simple  shear  tests  or  torsional  shear  tests  conducted 

under  cyclic  loading  conditions.  The  strain  amplitudes  in 

-2 

these  tests  may  vary  from  about  10  to  5  percent. 
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b)  forced  vibration  tests 

Forced  vibration  tests  may  be  performed  on  cylindrical  samples 

by  subjecting  them  to  longitudinal  vibrations  and  torsional 

#  2 

vibrations  with  strain  amplitudes  from  about  10  to  10 
percent.  The  tests  involve  the  determination  of  the  resonant 
frequencies  and  measurement  of  the  response  at  other  fre¬ 
quencies. 

c)  Free  vibration  tests 

In  these  tests  cylindrical  soil  samples  are  set  into  longitu¬ 
dinal  or  torsional  vibrations  as  in  the  forced  vibration  tests 
and  thereafter  the  decay  of  the  amplitudes  is  observed. 

These  tests  are  particularly  suited  for  the  determination  of 
the  damping  factor  since  the  logarithmic  decrement  at  any 
given  strain  level  is  obtained  directly  from  successive 
amplitudes  of  free  vibration. 

d)  Field  measurement  of  wave  velocities 

Field  tests  have  been  used  to  measure  the  velocity  of 

propagation  of  P-waves,  S-waves  and  Rayleigh  waves.  The  soil 

-4 

moduli  for  low  strains  of  about  5  x  10  percent  can  be 
determined  from  the  w' /e  . eloci^ies,  but  damping  factors  have 
not  been  obtained  from  these  tests. 


Linearity 


The  strain  amplitudes  in  soils  beneath  foundations  for  vibratory 
machines  are  usually  small  [42],  Test  data  for  soils  under  cyclic 
loading  conditions  have  shown  that  the  soil  response  is  approximately 

_3 

linear  for  small  strains,  say  10  percent  [43].  Thus  the  secant 
moduli  and  the  damping  factors  used  in  analyses  of  foundation  vibration 
problems  may  well  be  assumed  to  be  independent  of  the  strain  amplitude. 
Even  in  the  case  of  large  strains  which  are  developed  in  soils  during 
earthquakes,  researchers  have  successfully  employed  linear  analyses 
using  soil  properties  based  on  observed  strain  amplitudes  [20]. 
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2.4  Complex  Modulus  Representation 


Though  the  internal  damping  in  soils  is  not  considered  to  be 
caused  by  viscosity,  the  stress-strain  behavior  of  soils  under  vibratory 
loading  is  similar  to  that  of  viscoelastic  materials.  This  permits 
the  use  of  complex  moduli  in  representing  the  stress-strain  behavior 
of  soils  that  are  subjected  to  loads  varying  harmonically  in  time. 

The  stress-strain  relation  under  uniaxial  conditions  can  be  expressed  by 


in  which 


ac  exp  (iwt)  =  Ec  ec  exp  (icot) 


c  =  a  +  i  a2 


e  =  e1  +  i  c2 

EC(o)>  ^  E1(u)  +  i  E2  (w) 


(2.1) 

(2.2) 

(2.3) 

(2.4) 


are  the  complex  stress,  the  complex  strain  and  the  complex  modulus, 
respectively,  tu  is  the  circular  frequency  and  i  =  /*!.  The  superscript 

c 

denotes  a  complex  quantity. 

The  complex  stress  and  strain  may  be  visualized  as  a  pair  of 
vectors  rotating  at  the  frequency  w  about  the  origin  of  the  complex 
plane  *s  shown  In  fig.  8.  At  any  given  time,  the  actual  stress  and 
strain  are  the  projections  of  the  vectors  onto  the  real  axis.  The 
complex  modulus  is  generally  frequency  dependent.  Its  real  component, 
E^,  is  the  modulus  of  strain  which  is  in  phase  with  the  stress  and  its 
imaginary  component,  E2>  is  the  modulus  of  strain  which  is  90°  out  of 
phase  with  the  stress.  E^  is  associated  with  the  elastic  phenomenon  in 
which  energy  is  stored  in  a  recoverable  form  while  E2  is  associated 
with  the  viscous  phenomenon  in  which  energy  is  dissipated. 

Due  to  the  viscous  effect,  the  strain  vector  lags  behind  the 
stress  vector  by  an  angle  <f>T  which  lies  between  0°  and  90°.  This 

u 

angle  is  called  the  "loss  angle"  and  its  tangent  is  called  the  "loss 
tangent"  because  it  is  associated  with  energy  loss.  The  loss  tangent 
is  defined  by 

tan$L  =  E2/E1  (2.5) 

and  is  related  to  the  damping  factor  or  fraction  of  critical  damping. 
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6,  and  the  logarithic  decrement,  D,  by  [42] 


tand>.  =  2  8  ~  D/3.14 

Li 


(2.6) 


if  the  damping  is  small.  This  is  the  case  in  foundation  vibration 
problems. 

Since  most  data  on  damping  characteristics  of  soils  has  been 
presented  in  terms  of  the  fraction  of  critical  damping,  8  will  also  be 
used  here  as  the  measure  for  material  damping.  The  complex  modulus  may 
be  written  accordingly 


V  =  Ej  (1  +  i  2  8) 


(2.7) 


The  use  of  a  complex  modulus  implies  that  the  actual  hysteresis 
.loops  of  the  type  shown  in  Fig.  7  are  approximated  by  elliptical 
hysteresis  loops  which  are  equivalent  with  respect  to  the  slope  of  the 
principle  axes  and  the  areas  enclosed  by  the  loops.  Thus  the  secant 
modulus,  determined  as  shown  in  Fig.  7,  equals  the  real  component  of 
the  corresponding  complex  modulus,  and  the  energy  losses  per  cycle, 
which  are  proportional  to  the  areas  of  the  hysteresis  loops,  are  also 
matched. 

The  complex  modulus  EC  used  above  may  represent  any  of  the  usual 
moduli.  Young’s  modulus,  shear  modulus  or  bulk  modulus.  The  relations 
between  the  complex  moduli  of  an  isotropic,  linearly  viscoelastic 
material  at  a  given  frequency  are  analogous  to  those  between  the  real 
moduli  of  an  isotropic,  linearly  elastic  material.  Thus,  a  complex 
Lame's  constant  Xc  and  a  complex  Poisson's  ratio  vc  may  be  defined  by 


X  =  \1  +  i  X2  = 


GC  (EC  -  2  Gc) 
3GC  -  EC 


(2.8) 


v  =  v1  +  i  v2  = 


2  G 


-  1 


(2.9) 


in  which  EC  =  E^  +  i  E^  and  GC  =  +  i  Gj  denote  the  complex  Young's 
modulus  and  shear  modulus,  respectively. 

Values  for  Ec  and  Gc  can  be  directly  determined  from  longitudinal 
and  torsional  vibrations  of  cylindrical  soil  samples  in  forced  and  free 
vibration  tests.  However,  later  in  the  analyses  it  is  more  convenient 
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Co  u~e  the  complex  Lame’s  constants  \  and  G  • 

c 

Henceforth,  the  superscript  will  be  omitted  to  simplify  the 
notation.  It  is  understood  that  X  and  G  denote  complex  moduli  unless 
otherwise  stated. 


2. 


Wave  Velocities 


The  velocity  of  P-waves,  Vp,  and  the  velocity  of  S-waves,  Vg,  in 
linearly  viscoelastic  media  are  [17] 


Vp  =  |Re  /p/(X  +  2G) 

V$  =  (Re  v^Tg]"1 

l  1 


-1 


(2.10) 

(2.11) 


in  which  0  is  the  mass  density  and  the  symbol  Re  specifies  the  real 
part.  Vp  and  Vg  vary  with  frequency  as  do  the  complex  Lame's  consts 
If  the  imaginary  parts  of  the  complex  Lame's  constants  are  small. 


Vp  =  Re  /(X+2G)/p  and  Vg  =  Re  /G/p. 

In  linearlv  elastic  media,  the  P-wave  and  S-wave  velocities  are 
independent  of  frequency.  They  are 

(2.12) 


Vp  =  /(X  +  2G)/p 


\j  —  Jn  /*-. 
•s  _  tG,v 


in  which  G  and  X  denote  real  Lame's  constants. 


/  O  1  *>  N 

JL.); 


2.6  Material  Damping  in  Foundation  Vibrations 

The  frequency  of  motion  hardly  affects  the  secant  moduli  and  the 
damping  values  of  soils  throughout  the  practical  frequency  range. 
Therefore  the  complex  moduli  may  be  considered  frequency  independent. 
This  means  that  the  damping  forces  in  soils  are  proportional  to  the 
strain  amplitudes  rather  than  to  the  strain  velocities  [43]. 

The  fraction  of  critical  damping  in  foundation  vibration  problems 
is  usually  small.  Its  values  for  shear  deformations  may  typically 
vary  from  2  to  5  percent,  and  its  values  for  compressi onal  deformations 
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3.  METHOD  OF  ANALYSIS 


3. 1  Method  of  Complex  Variables  for  Harmonic  Motion 


It  is  convenient  to  represent  the  forces,  displacements,  strains 
and  stresses  of  harmonic  motion  by  complex  variables.  If  f(t)  is 
harmonic  in  time  it  can  be  defined  by 


f(t)  =  f  cos(wt  +  (j)) 

(3.1) 

or 

f(t)  =  f^  coscot-f2  sinut 

(3.2) 

or 

f(t)  =  Re  (fC  exploit) 

(3.3) 

if.  which  o>  is 

c 

the  circular  frequency  and  f  is 

the  complex  variable 

fC  =  fx  +  i  f2 

(3.4) 

with  i  =  /-I. 

The  amplitude  f  is 

£o  =  V£l+f2 

(3.5) 

and  the  phase 

angle  <j>  is 

4>  =  tan'1(f2/f1) 

(3.6) 

Equation  of  Motion  for  Discrete  Systems 

The  equation  of  motion  for  a  simple  damped  oscillator  with  a  mass 
M,  a  spring  constant  K,  and  a  dashpct  constant  C  is 

K  u(t)  +  C  u(t)  +  M  u(t)  =  P(t)  (3.7) 

m  which  u,  u,  and  u  are  the  displacement  velocity  and  acceleration  of 
the  mass,  respectively,  and  P(t)  is  the  driving  force.  For  harmonic 
motion  at  the  frequency  w,  Eq.  (3.7)  can  be  reduced  to 

(Kc-w2  M)  uc  =  Pc  (3.8) 
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in  which 

KC=K1+iK2=K+iCu  (3.9) 

uC  -  u^  +  i  u^  (3. 10) 

PC  =  Px  +  i  P2 

Eq.  (3.8)  is  a  time  independent,  complex,  linear  equation  which  is 
equivalent  to  the  two  real  equations 

Kx  u1  -  K2  u2  -  u)2  M  u1  =  P  (3.12) 

^1  u2  +  ^2  U1  ”  ^  ^  u2  =  ^2  (3.13) 

The  first  equation  states  equilibrium  when  sinu)<:=0,  the  second  when 
cosut=0.  If  both  equations  are  satisfied  then  equilibrium  prevails  at 
any  time. 

The  equations  of  harmonic  motion  for  a  system  of  n  degrees  of 
freedom  are,  in  generalization  of  Eq.  (3.8), 

([Kc]  -  a)2  [M])  {uc}  =  {Pc}  (3.14) 

in  which  the  mass  matrix  [M]  and  the  complex  stiffness  matrix  [Kc]  =  [K^] 
+  i[K2]  are  of  order  n  x  n.  [K^]  contains  the  stiffness  coefficients 

and  [K2]  the  damping  coefficients  multiplied  by  w.  The  vectors  (uC)  and 

f  C  i  c  c 

IP  )  contain  the  complex  displacements  and  forces,  u.  and  P.,  j=l,  ...» 

c  3  ^ 

n,  respectively.  The  complex  stiffness  matrix  [K  ]  may  be  generated  in 
the  same  manner  as  a  real  stiffness  matrix,  the  only  difference  being 
that  complex  moduli  are  used  instead  of  real  moduli. 


Equation  of  Motion  for  Continuous  System 

The  equilibrium  conditions  for  an  infinitesimal  element  together 
with  the  strain-displacement  relations  and  the  stress-strain  relations 
yield  the  field  equations  of  motions  in  terms  of  the  displacements. 

If  the  motion  is  harmonic  the  field  equations  can  be  reduced  to 
time  independent  equations  in  the  same  way  as  shown  above.  Using 
rectangular  Cartesian  coordinates,  x,  y,  z,  as  indicated  in  Fig.  4,  and 
assuming  plane  strain  conditions,  i.e.  all  derivatives  with  respect  to 
y  vanish,  the  field  equations  for  the  homogeneous  case,  in  which  no 
external  body  forces  act  on  the  element,  are 
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imstm 
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4-  (X  +  G) 


+  PU)2  u£  =  0 


(3.15a) 


/32uC 

ll  1 


^3x 


^2  C\  /« 2  c  cx 

3  uz\  6  ux  3  M 

TF)  +  <X  +  G>  VS*  +  ^f) 


c 

o  u 


v3x 


~2 
3  u 


3z 


+  pu2  uC  =  0 
z 


,  2  c  „ 

+  pw  u  =0 

y 


(3.15b) 


(3.15c) 


C  (.  r 

u^,  u.  ,  and  are  the  complex  displacements  in  the  x»  y,  and  z- 
dirc.tions,  respectively,  p  is  the  mass  density  and  X  and  G  are  the 
complex  moduli  introduced  in  Section  2.2.  Eqs.  (3.15)  are  valid  for 
small  deformations  in  an  isotropic,  linearly  viscoelastic  medium.  Eqs. 
(3.15a)  and  (3.15b)  are  coupled.  They  govern  the  motion  in  the  x-z- 
plane,  while  Eq.  (3.15c)  governs  the  motion  in  the  y-direction. 

The  analogous  equations  for  axisymmetric  conditions,  when  all 
derivatives  with  respect  to  the  angular  direction  9  vanish,  are 


(X  +  2G) 


/32uC  ,  3uC  uC  3uC\  / 32uC  32uC\  2  c  _n 

{ _ r  1 _ r _ r  .  _ z  }  ,  r  ( _ r _ z  )  +  pw  u  =  0 

\3r2  r  3r  r2  3r3z/  '3z2  3r3z' 


(3.16a) 


(^2  c  -  c  ,2  c\  /-.2  c  -  c  ~2  c  c> 

o  u  ,  du  d  u  \  /3  u  ,  du  d  u  .  du 

^  +  7  3T  ^j+ G  +  7  3T  -  3^  -  7  ^ 


+  pLJ2  UC  =  0 
r 


(3.16b) 


/^+iK 

\^2  r3r  2  -.2/ 

'Sr  r  3z  > 


c 

'e 


=  o 


(3.16c) 


m  which  u^,  Ug,  and  are  the  complex  displacements  in  the  r,  9,  and 
2-directions  of  a  cylindrical  coordinate  system  as  indicated  in  Fig.  5. 

Eqs.  (3.16a)  and  (3.16b)  are  coupled.  They  govern  the  motion  in  the 
r-z-plane,  whereas  the  motion  in  the  6-direction  is  governed  by  Eq.  (3.16c) 
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It  is  expedient  to  have  a  stationary  principle  available  for  the 
viscoelastic  medium  as  discretized  by  the  finite  element  method.  To 
this  end,  the  well  known  principle  cf  virtual  work  will  be  specialized 
for  application  to  harmonic  motion.  The  principle  of  virtual  work  for 
the  dynamic  case  can  be  stated  in  the  following  form  [17] 


f  {c}T{o)dV  -  f  (u}T({f)  -  p(u})dV  -  f  {u}T{T)dS  =  0  (3.17) 

V  v  JS 

in  which  the  stress  vector  {o},  the  body  force  vector  (f),  the  surface 
force  vector  {t},  and  the  acceleration  vector  {ii}  are  functions  of  space 
and  time.  The  vectors  {u}  and  (e }  contain  the  virtual  displacements  and 
virtual  strains,  respectively.  The  integration  /^  . . .  dV  is  performed 
over  the  total  volume  of  the  body  and  the  integration  ...  dS  is  per¬ 
formed  over  that  part  of  the  surface  of  the  body  where  no  displacements 
are  prescribed.  The  virtual  displacements  are  small,  triply  differen¬ 
tiable  functions  of  space  and  must  vanish  over  the  surface  of  the  body 
where  displacements  are  prescribed  [17].  Accordingly,  the  virtual 
strains  are  small,  twice  differentiable  functions  of  space. 

Eq.  (3.17)  states  that  equilibrium  prevails  if  the  sum  of  the 
virtual  work,  performed  by  the  actual  forces  and  actual  strains  on  the 
virtual  displacements  and  virtual  strains,  vanishes. 

Alternatively,  the  rate  of  the  virtual  work  may  be  considered. 

This  leads  to  the  equation  of  motion  in  the  following  stationary  form 


f  {e)T{o}dV  -  /  {u}T({f}  -  p{ii})dV  -  /  {{i}T{T}dS  =  0  (3.18) 

V  '  V  JS 

A  A 

in  which  the  vectors  {u}  and  {e}  contain  the  virtual  velocities  of  the 
displacements  and  strains,  respectively. 

If  the  motion  is  harmonic  with  the  circular  frequency  w,  each 
variable  in  Eqs.  (3.17)  and  (3.18),  including  the  virtual  displacements 
and  the  virtual  strains,  can  be  described  in  time  as  is  f(t)  in  Eq. 
(3.2).  The  time  derivatives  in  Eqs.  (3.17)  and  (3.18)  are 
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(3.19) 


{ li }  =  -  co  ({u^}  costot  -  {t^}  sintot) 


{ u }  =  -  to  ({u^}  sincot  +  (u2J  costot) 

A 

{£}  =  -  to  ({e^}  sincot  +  {z^}  costot) 


(3.20) 


(3.21) 


The  integration  of  Eqs.  (3.17)  and  (3.18)  over  the  time  of  one  period 
leads  to 


/  [{ei}T{a1)  +  {e2)T(o2}  -  {ai}T({F1>  +  pto2{Ul»  -  {u2)T((F2}  +  pco2{u2})]dV 


-  f  ({G1>t{t1  }  +  (u2HT2})dS  =  0 


(3.22) 


/  [{e1)T{a2)  -  {e2}T{Pl}  -  {G1)T({F2>  +  Pto2(u2))  +  {G2}T({f1>  +  Pto2{u1»]dV 
-  /  ({uj1'^}  -  {u2}T{T1»dS  =  0  (3.23) 


respectively,  since 


/tt/io  2  /-27r/co  2 

sin  tot  dt  =  J  cos  cot  dt  =  tt/co 


--21T /to 

/  sincot  •  coscot  dt  =  0 


(3.24) 


(3.25) 


Now,  Eq.  (3.23)  is  multiplied  by  i  =  /-I  and  added  to  Eq. 
(3.22).  If  complex  variables  are  introduced  according  to  Eqs.  (3.3)  and 
(3.4),  this  leads  to  the  equation  of  harmonic  motion  in  the  stationary 


f  {€>  {aC)dV  -  /  {G}  ({FC}  +  Pco2  (u  })dV  -  /  {G}  (TC}dS  =  0 

Jy  V  S 

(3.26) 

£ 

in  which  the  superscript  indicates  a  complex  variable,  as  before,  and 
the  indicates  the  transposed  complex  conjugate  (e.g.  (uC}={u^}+i{u2}and 


{G}*={G1}T-i{u2}T). 
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Henceforth,  however,  the  superscript  c  will  be  omitted  for  con¬ 
venience,  and  it  is  understood  that  the  quantities  of  harmonic  motion 
are  represented  by  complex  variables  as  is  f(t)  in  Eq.  (3.3). 


3.2  Finite  Element  Method 


The  finite  element  method  is  a  numerical  technique  for  obtaining 
approximate  solutions  to  complex  boundary  value  problems.  The  develop¬ 
ment  of  the  method  in  the  field  of  structural  mechanics  began  in  the 
early  1950's  with  the  work  of  Turner  et  al.  [45].  Since  then  the  method 
has  been  developed  extensively  and  is  now  widely  used  in  structural 
and  continuum  mechanics.  It  has  also  been  applied  successfully  to 
various  other  physical  problems. 

A  recently  published  text  by  Zienkiewicz  [53]  describes  in  detail 
the  application  of  the  finite  element  method  in  structural  and  continuum 
mechanics.  The  text  also  contains  an  extensive  list  of  references  on 
the  subject.  The  finite  element  method  has  been  used  very  successfully 
in  static  and  dynamic  analyses  of  plane  and  axisymmetric  problems  [12, 

13,  50,  51].  It  has  also  been  applied  to  quasi-static  and  dynamic 
analyses  of  viscoelastic  continua  with  considerable  success  [1,  3,  11, 
32,  35]. 

Since  the  finite  element  method  is  widely  known,  only  the  principal 
steps  involved  in  the  particular  application  of  the  method  to  the 
analysis  of  harmonic  motion  in  a  viscoelastic  medium  are  summarized 
below. 


(i)  Discretization: 


The  irregular  region  I  is  subdivided  into  subregions,  called 
finite  elements,  as  shown  in  Fig.  6.  The  finite  elements  are  prisms  in 
the  plane  case  and  rings  in  the  axisymmetric  case  and  have  arbitrary 
quadrilateral  cross-sections.  They  are  interconnected  at  straight  or 
circular  nodal  joints,  which  are  located  at  the  edges  of  the  elements. 
Because  the  strains  are  either  plane  or  axisymmetric,  the  displacements 
vary  only  over  the  cross-section  of  the  elements.  It  is  assumed  that 
the  displacement  field  of  an  element  is  restricted  to  a  finite  number 
of  degrees  of  freedom.  Here,  the  displacements  of  the  nodal  joints  are 
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introduced  as  the  degrees  of  freedom.  Each  degree  of  freedom  of  an 
element  is  associated  with  a  shape  function  which  defines  the  displace¬ 
ment  variation  within  the  element.  The  shape  functions  vary  linearly 
along  the  boundary  of  an  element,  see  Section  3.3.  Hence,  the  displace¬ 
ments  across  the  boundary  between  adjacent  elemer.  *s  remain  continuous 
during  deformation  if  the  elements  are  interconnected  at  their  nod.al 
joints. 


(ii)  Element  matrices: 

The  influence  coefficients  of  an  individual  element  with  respect 
to  unit  nodal  displacements  are  derived  from  the  principle  of  virtual 
work,  Eq.  (3.26).  This  step  is  performed  in  Section  3.3.  The  influence 
coefficients  are  collected  in  the  element  stiffness  matrix  [K * ]  and  the 
element  mass  matrix  [M*].  Because  harmonic  motion  of  a  viscoelastic 
body  is  considered,  IK']  is  here  a  complex  matrix  representing  the 
viscous  as  well  as  the  elastic  resistance  of  the  element  against 
deformation.  The  matrix  [M*]  is  real  and  represents  the  inertia  forces 
of  the  element. 

(iii)  Assemblage  of  finite  elements: 

The  elements  are  assembled  by  matching  the  displacements  of 
adjacent  elements  at  common  nodal  joints.  The  global  stiffness  and  mass 
matrices  of  the  assemblage,  [K]  and  [M] ,  are  accordingly  formed  by 
addition  of  the  element  matrices  while  observing  the  global  numbering 
system  of  the  nodal  displacements.  The  procedure  is  that  of  the  direct 
stiff  r.ess  method  [53]. 


(iv)  Equation  of  motion: 

Nodal  displacements  are  prescribed  which  provide  the  kinematic 
stability  of  the  body  and  external  loads  are  applied  as  discrete  forces 
acting  at  the  nodal  joints  or  as  prescribed  nodal  displacements.  The 
loads  vary  harmonically  with  time  at  the  circular  frequency  co.  The 
equation  of  motion  is  therefore 


in  which 


[A]  {u}  =  {bl 


[A]  =  [K]  -  u>  [M] 


(3.27) 


(3.28) 
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is  a  symmetric  banded  matrix  with  complex  coefficients.  The  vector  {u} 
contains  the  complex  nodal  displacements  and  the  vector  {b}  the  complex 
nodal  forces. 

Eq.  (3.28)  can  be  partitioned  in  the  following  way 


Aff  Afs  Uf 

aE  A  u 

fs  ssJ  Is 


(3.29) 


in  which  (u^)  contains  the  unconstrained  or  free  nodal  displacements  and 
{ ug }  the  prescribed  nodal  displacements.  Accordingly,  {b^}  contains 
the  prescribed  nodal  forces  and  (bg)  the  reaction  forces.  Hence,  Eq. 
(3.29)  can  be  reduced  to 


lAff]  (uf >  =  (bf)  -  [Afs]  (Ug > 


(3.30) 


(v)  Solution  of  equation: 

Eq.  (3.30)  constitutes  a  set  of  simultaneous  complex  linear 
equations.  The  coefficient  matrix  is  symmetric,  sparse  and  also  narrowly 
banded  if  a  favorable  numbering  system  of  the  nodal  displacements  is 
used.  The  equations  may  be  solved  by  a  Gaussian  elimination  algorithm 
without  pivots  using  complex  arithmetic.  Once  the  displacement  vector 
[uf)  is  known  the  reaction  forces  {bg}  can  be  computed  by 


{bs}  =  lAfs]  {uf}  +  tAss3  {us} 


(3.31) 


Further  details  of  the  procedure  outlined  above  are  given  in 
Ref.  [53] .  However,  a  few  points  need  additional  discussion. 

The  semi-infinite  layered  regions  R  and  L  must  also  be  considered. 
Dynamic  stiffness  matrices  for  these  regions  are  derived  in  Chapters  4 
through  7  for  the  cases  of  plane  and  axisymmetric  Love  and  Rayleigh  wave 
motion.  When  these  matrices  are  added  by  the  direct  stiffness  method 
into  the  global  matrix  [A]  of  Eq.  (3.27),  the  influence  of  the  semi¬ 
infinite  regions  on  the  motion  in  the  irregular  region  is  properly 
accounted  for. 

The  accuracy  of  the  analysis  depends  upon  the  fineness  of  the 
finite  element  mesh.  In  general,  the  larger  the  strain  gradient  the 
finer  should  be  the  mesh.  Therefore  the  mesh  should  be  refined  in  areas 
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where  large  strain  gradients  are  to  be  expected. 

Characteristic  features  of  harmonic  motion  in  elastic  or  visco¬ 
elastic  media  are  the  wave  lengths  of  P  and  S-waves.  In  order  to  obtain 
a  good  approximation  to  the  actual  wave  motion,  the  maximum  cross- 
sectional  dimension  of  any  finite  element  should  be  small  compared  to 
the  length  of  S-waves,  which  is  always  shorter  than  the  length  of  P- 
waves  in  the  same  medium  at  any  given  frequency.  Lysmer  and  Kuhlemeyer 
[25]  recommended  that  the  maximum  dimension  of  a  rectangular  finite 
element  with  a  linear  strain  field  be  not  greater  than  about  1/10  of 
the  length  of  S-waves  in  the  medium.  This  recommendation  also  applies 
to  the  quadrilateral  finite  elements  presented  in  Section  3.3. 

Numerical  solutions  are  compared  with  analytical  solutions  in 
Chapters  9  and  11 .  The  comparisons  show  clearly  that  very  good  results 
can  be  obtained  by  the  proposed  method  for  steady-state  wave  propagation 
problems  in  semi-infinite  elastic  and  viscoelastic  media. 

It  was  mentioned  earlier  that  the  set  of  complex  linear  equations, 

Eq.  (3.30),  may  be  solved  by  a  Gaussian  elimination  algorithm  without 
pivots.  The  algorithms,  used  in  the  computer  programs  of  Appendix  5  and 
6,  are  coded  in  the  FORTRAN  IV  language  and  employ  the  complex  arithmetic 
capabilities  of  FORTRAN  IV.  By  utilizing  complex  arithmetic,  computer 
storage  requirements  are  half  that  required  for  the  equivalent  real 
arithmetic  computation,  in  which  the  complex  matrix  and  the  complex 
sectors  are  decomposed  into  real  and  imaginary  parts.  Since  the  co¬ 
efficient  matrix  is  banded  and  symmetric,  only  half  of  the  band  needs 
to  be.  stored  in  the  computer  memory. 

No  numerical  sensitivity  has  been  noted  in  the  solution  of  the 
equations.  Numerical  sensitivity  could  arise  when  the  frequency  of 
excitation  is  extremely  close  to  the  natural  frequency  of  a  subsystem 
of  the  discretized  model.  But  this  cannot  be  the  case  if  the  model 
contains  some  viscous  damping  associated  with  each  degree  of  freedom, 
because  then  none  of  the  subsystems  possesses  a  real  natural  frequency 
whereas  the  frequency  of  excitation  is  always  real.  If  the  model  contains 
no  damping,  the  natural  frequency  of  a  subsystem  could  coincide  with 
the  exciting  frequency.  But  this  is  very  unlikely.  Therefore  an  al¬ 
gorithm  without  pivots  is  used  and  the  extra  computation  time  and  storage 
needed  for  pivoting  are  saved. 
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3.3  Element  Stiffness  and  Mass  Matrices 


In  order  to  assemble  the  matrices  [K]  and  JM]  appearing  in  Eq.  (3.28) 
it  is  necessary  to  first  determine  the  matrices  IK* 3  and  [M')  for  the 
individual  elements  of  the  assemblage.  The  plane  or  axisymmetric  cross- 
section  of  a  typical  element  is  shown  in  Fig.  9,  which  also  shows  the 
global  coordinate  system  (r,  z)  and  a  local  coordinate  system  (n,  O • 

The  two  sets  of  coordinates  are  related  through  the  transformation 


<r  z>  =  {b} 


(3.32) 


in  which  r  ,  z  ,  m=l,2,3,4,  are  the  coordinates  of  the  four  nodal  points 
mm 

of  the  element  cross-section  and  {b}  is  the  column  vector  defined  by 


(l-o  (1-n) 

1  d-o  (i+o) 

{b}  =  4  (1+0  (1+0) 


(3.33) 


[d+O  (i-n)  J 

The  Jacobian  matrix  of  the  transformation  is 

fri  zil 


[J]  - 


=  ID] 


r2  Z2 


r3  z3 


r4  Z4 


(3.34) 


in  which 


[D]  =  T 


-(i-n)  -d+n)  (l+n)  (l-n) 


i  r-(i- 

«  l-a- 


o  (l-o  (l+o  -(1+0 


(3.35) 


The  determinant  of  [J]  is  denoted  J. 

For  both  the  Love  wave  case  and  the  Rayleigh  wave  case  the 
assumption  will  be  made  that  the  displacements  within  the  element  vary 
according  to 


<5(r,z)  =  S(£,n) 


(3.36) 


(b}T 


in  which  5^,  6^,  5^  and  6^  are  the  nodal  displacements  in  some 
coordinate  direction  of  the  element  shown  in  Fig.  10.  Due  to  the  special 
properties  of  the  chosen  transformation,  Eq.  (3.32),  this  implies  that  all 
displacements  vary  linearly  along  the  boundaries  of  the  element  and  that 
all  elements  remain  compatible  as  the  nodal  joints  are  displaced. 


Element  Matrices  for  Axisymmetric  Love  Wave  Case 

In  the  Love  wave  case,  the  quadrilateral  element  shown  in  Fig.  10a 
has  four  degrees  of  freedom  which  arc  the  nodal  displacements  in  the 
O-direction.  The  displacement  field  6g  =  o(£,  r|)  is  defined  by  Eq. 
(3.36). 

The  strain-displacement  relations  in  cylindrical  coordinates  are 


96  £ 

Yr0  9r  “  r 

Differentiation  by  the 


,  36 

and  *,e  ■  3? 

chain  rule  yields 


36 

3r 

<  ► 
36 
3z 


[T]  {6} 


in  which  [T]  is  the  2x4  matrix 


(3.37) 


(3.38) 


[T]  =  [J3_1[D]  (3.39) 

and  the  matrices  [J]  and  [D]  are  those  defined  in  Eqs.(3.34)  and  (3.35). 
The  strains  expressed  in  terms  of  the  nodal  displacements  are 


(e) 


=  [S]  {6} 


(3.40) 
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in  which  [S]  is  the  2x4  matrix 


[s]  = 


'll  -  Vr 


b12  -  b2/r 


'l3  '  b3/r 


b14  ‘  b4/r 


(3.41) 


L  C21  t22  C23  C24  J 

The  elements  t  ,  i*l,  2,  j=l,  2,  3,  4,  are  the  elements  of  [T]  and 
j=l,  2,  3,  4,  are  the  elements  of  (b),  defined  by  Eq.  (3.33). 
The  stress-strain  relation  of  an  isotropic  material  is 


{a} 


=  G{e} 


(3.42) 


in  which  G  represents  the  complex  shear  modulus  if  harmonic  motion  in  a 
viscoelastic  material  is  considered. 

The  nodal  forces  Q.,  j=l,  2,  3,  4,  which  are  shown  in  Fig.  10a  and 
are  collected  in  the  vector  (q),  are  the  only  external  forces  acting  on 
the  element.  They  vary  harmonically  in  time.  Since  [S]  and  {b}  are 
real,  the  principle  of  virtual  work,  Eq.  (3.26),  yields 

G  f  (6}*[S]T[S]{6}dV  -  obi2  f  {6}*{bHb}T{6}dV  »  {6}*{Q}  (3.43) 

V  V 

A 

in  which  the  vector  {6}  contains  the  virtual  displacements  of  the  nodal 
* 

joints  and  the  denotes  the  conjugate  transposed. 

The  four  virtual  displacement  states  in  each  of  which  a  single 

A 

component  of  {6}  is  assigned  the  value  1  while  the  other  three  com¬ 
ponents  are  0  yield  the  matrix  equation 

([K’J  -  u)2[M']){6}  =  (q)  (3.44) 

in  which  [K’]  is  the  4x4  element  sttrfness  matrix 

Ik* ]  =  G  Jv  [S]T[S]dV  (3.45) 

and  [M’3  is  the  4x4  consistent  element  mass  matrix  15] 

[M’]  =  p  fv  (bHb)TdV  (3.46) 


The  volume  integral  can  be  converted  to  an  area  integral  since 
dV=r  d8  dA  for  an  axisymmetric  element.  In  the  £-r,  coordinate  system 
dA=J  d£  dn*  The  stiffness  and  mass  matrices  for  a  one  radian  segment 
of  the  element  are  therefore 
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mamsmi 


(3.47) 


[K* ]  =  f  [  G[S]T[S]  r  J  dC  dn 

-1  -1 

[M* ]  =  f  f  p{b}{b}T  r  J  d(  dn  (3.48) 

-1  -1 

The  integration  is  performed  numerically  by  Gaussian  quadrature 
with  a  four  point  scheme  [53] . 

Element  Matrices  for  Plane  Love  Wave  Case 

Axisymmetric  strain  conditions  become  plane  strain  conditions  at 
large  distance  from  the  axis  of  symmetry.  The  matrices  for  a  plane 
strain  element  can  therefore  be  deduced  from  Eqs.  (3.47)  and  (3.48)  by 
considering  a  segment  of  infinite  radius  subtending  a  unit  arc. 

Since  [3]  reduces  to  [T] ,  defined  by  Eq.  3.39,  the  matrices  are 


[K']  =  f1 

f  g[t]t[t]  j  d£  dn 

(3.49) 

J-1 

-l 

f  1 

r 1  T 

[M’J  =  I 

/  p{bHb}  J  dt  dn 

(3.50) 

-1 

4 

Element  Matrices  for  Axisymmetric  Rayleigh  Wave  Case 


The  Rayleigh  wave  case  involves  two  displacements  and  two  forces 
per  nodal  joint  as  shown  in  Fig.  10b. 

The  r  and  z  displacements  within  the  element  are,  according  to 
Eq.  (3.36), 


r5r 

1/  =  IN]{6} 

l  z 


(3.51) 


in  which  {6}  contains  the  nodal  displacements  6^,  j=l,  ...,  8,  and  [N] 
is  the  2x8  matrix 
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[N]  = 


bL  0  b,  0  b3  0  b4  0 

0  b.  0  b0  0  b,  0  b, 
12  3  4 


(3.52) 


The  elements  b^ ,  j=l*  2,  3,  4,  are  the  components  of  \b}  defined  in 
Eq.  (3.33). 

The  strain-displacement  relations  are 


;  e  =  3— 
z  3z 


(3.53a) 


6  36  36 

r  r  ,  z 

£0  r  ;  Yrz  "  3z  +  3r 


(3.53b) 


Differentiation  by  the  chain  rule  leads  to 


(e)  = 


=  [S] {6} 


(3.54) 


in  which  [$'•]  is  now  the  4x8  matrix 


C11  0  t12  0  C 13  °  t14  0 


IS]  = 


0  C21  0  C22  ° 


t23  ° 


b  /r  0  b„/r  0  b_/r  0  b./r  0 
12  3  4 


(3.55) 


[.  t21  tll  C22  t12  t23  t13  fc24  fc14  J 

The  elements  t„,  i=l>  2,  j=l,  2,  3,  4,  are  the  elements  of  [T]  defined 
in  Eq.  (3.39)  and  b  ^ ,  j=l,  2,  3.  4,  are  the  elements  of  {b}  defined  in 
Eq.  (3.33). 

The  stress-strain  relation  for  an  isotropic  material  is 


(a)  = 


=  [C]  {e } 


(3.56) 


iflliM  Wi ,  t :  u  i>  finim  Mint'll, 'if  ‘I  liwll/Mi'iJ-tiUii1  fil , 


in  which 


X  +  2G  X  X  0 

X  X  +  2G  X  0 
ICl  X  X  +  2G  0 

0  0  0  G 

If  harmonic  motion  in  a  viscoelastic  material  is  considered,  the  Lame’s 
constants  X  and  G  are  complex  moduli,  which  are  generally  frequency 
dependent. 

The  further  derivation  of  the  stiffness  and  mass  matrices  is 
analogous  to  that  in  the  Love  wave  case  and  leads  to  the  8x8  matrices 
for  a  one  radian  segment  of  the  element 


IK’]  =  f1  fl[S]T[C][S]  r  J  dg  dn 

(3.58) 

-1  J-l 

,  1  ,  1 

[M’l  =  /  /  olN]  [N]  r  J  dC  dn 

(3.59) 

-1  -1 

in  which  [S]  is  the  matrix  defined  by  Eq.  (3.55). 

The  integration  is  again  performed  numerically  by  Gaussian 
quadrature  with  a  four  point  scheme. 


Element  Matrices  for  Plane  Rayleigh  Wave  Case 

The  stiffness  and  mass  matrices  for  the  plane  strain  element  can 
be  deduced  from  those  for  the  axisymmetric  element  by  considering  an 
element  segment  of  infinite  radius  subtending  a  unit  arc.  This  leads 


IK’]  =  rVlCHSJ  JdCdn 

-l  -l 

r  1  r 1  t 

[M']  =  /  /  PtNriNj  J  d£  dn 

-1  -1 


(3.60) 


(3.61) 


in  which  [3]  is  the  3x8  matrix  obtained  by  deleting  the  third  row  in 
matrix  (S]  and  [C]  is  the  3x3  matrix  obtained  by  deleting  the  third 
row  and  the  third  column  in  raacrix  ICj . 
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4.  PLANE  LOVE  WAVE  MOTION 


4 . 1  Eigenvalue  Problem  for  Continuous  Layered  Region 

Free  harmonic  motion  under  plane  strain  conditions  in  a  semi¬ 
infinite  layered  region  of  the  type  shown  in  Fig.  4  is  considered.  The 
region  consists  of  n  isotropic  linearly  elastic  or  viscoelastic  layers, 
which  are  welded  together  at  their  interfaces.  The  jth  layer  has  the 
mass  density  p..  and  the  complex  shear  modulus  .  All  displacements  are 
perpendicular  to  the  x-z-plane  and  are  described  by 

6  =  6^(x,z,t)  =  u(x,z)  exp(iu)t)  (4.1) 

in  v* ich  w  is  the  circular  frequency. 

The  governing  homogeneous  differential  equations  for  the  spatial 
part  of  the  displacements  in  the  n  layers  are,  according  to  Eq.  (3.15c), 


+  Pj  to  u  =  0 


j=l> 


(4.2) 


These  partial  differential  equations  may  be  reduced  to  ordinary 
differential  equations  by  separation  of  the  variables.  The  assumption 
that 

u(x,z)  =  v(z)  *  g(x)  (4.3) 

leads  to 


2-1  A  +  °-± 

.  2  v  G. 

VZ  J 


0) 


2  g 


3x 


j=l,  n 


(4.4) 


which  can  only  be  valid  for  all  values  of  x  and  z  in  the  respective 
layer  if  both  sides  of  Eq.  (4.4)  are  equal  to  the  same  constant.  Setting 


this  constant  equal  to  k  yields  the  ordinary  differential  equations 

j=l,  ...»  n  (4.5) 


2 

+  (<c2p./G.  -  k2)  v  =  0 
dz  J  J 
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MS 


and 


^f+k2g  -  0 

dx 


(A.  6) 


The  latter  equation  has  a  solution 

g  =  exp(-ikx)  (4.7) 

in  which  i  =  /-I.  The  parameter  k  is  called  the  wave  number,  which,  if 
real,  is  related  to  the  phase  velocity,  c,  and  to  the  apparent  wave¬ 
length  in  the  x-direction,  1,  by 

(4.8) 


,  OJ  27T 

k"  c*T 


Thus,  a  solution  to  Eq.  (4.2)  may  be  expressed  in  the  form  of  a  wave 

6  =  v(z)  exp(iwt-ikx)  (4.9) 

The  function  in  the  z-direction,  v(z) ,  which  will  be  called  the  mode 
shape  of  the  wave,  must  satisfy  the  n  ordinary  differential  equations, 
Eq.  (4.5).  In  addition,  the  boundary  conditions  at  the  planes,  z=zy 
j=l,  ...,  n+1,  i.e.  continuity  of  the  displacements  and  shear  stresses 
at  the  interfaces  between  the  n  layers,  zero  displacements  at  the  rough 
rigid  base  and  zero  shear  stresses  at  the  free  surface,  must  also  be 
satisfied.  Since  the  shear  stress  on  a  horizontal  plane  is 


t  =  G  -r—  =  G  ~  exp(iwt-ikx) 
zv  dz  dz  r 


these  boundary'  conditions  may  be  expressed  in  the  form 

j— 2, 3 ,  . .  * ,  n 


v(z!)  =  v(zV) 
J  1 


G.  4~  v(z’)  =  G.  ,  4~  v(zV) 
3  dz  3  3+1  dz  3 


j=2,3,  ...n 


(4.10) 

(4.11a) 

(4.11b) 


and 


tz  V(z=0)  =  0 


v(z-H)  =  0 


(4.11c) 


(4. lid) 


in  which  zj  and  z?  refer  respectively  to  planes  immediately  above  and 


immediately  below  the  interface  z y 
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The  general  solution  to  Eq.  (4.5)  for  z^<z<z^+1  is 

v(z)  =  a  cosy.z  +  b.  sinY.z 
3  3  3  3 


in  which 


7  7  7 

YV  =  w  P-/G.  -  k 

3  3  3 


(4.12) 


(4.13) 


The  shear  stresses  at  the  interfaces  j  and  i+1  in  terms  of  v.=v(z.)  and 

3  3 

Vj+j=v(zj+^)  follow  from  Eqs.  (4.10)  and  (4.12)  and,  after  some  algebra, 
reduce  to 


-T  (zV)  =  (c  v  +  d  v  )  exp(iwt-ikx) 
zy  J  33  3  3+1 


and 


Tzy(zj+1)  =  (djvj  +  CjVj+l)  ^(^t-i-lcx) 


in  which 


and 


c.  =  +  G.y.  cotYjhj 

d.  =  -  G.Y./siny.hj 


h.  =  z. ,,-z. . 
3  3+1  3 


(4.14a) 

(4.14b) 

(4.15a) 

(4.15b) 


Substitution  of  Eqs.  (4.14)  into  Eqs.  (4.11)  leads  to  a  set  of  n 
homogeneous  transcendental  equations 


C1V1  +  d2V2 


=  0 


d.  -v.  +  (c.  .  +  c.)  v.  +  d.v  =  0 

3-1  3-1  3-1  3  3  3  3+1 


3  *  *  *  j  u 


d  ,v  .  +  (c  -  +  c  )  v 
n-1  n-1  n-1  n  n 


=  0 


(4.16) 


Non-trivial  solutions  to  these  equations  can  be  obtained  only  for  eigen- 
2 

values,  k  ,  which  cause  the  determinant  of  the  coefficient  matrix  to 
vanish.  Eq.  (4.16)  states  a  difficult  eigenvalue  problem,  because  the 
eigenvalues  are  "hidden"  in  the  arguments  of  the  transcendental  ex¬ 
pressions  c.  and  d..  Infinitely  many  different  eigenvalues  and  associated 
eigenfunctions  exist;  but  it  is  difficult  to  obtain  even  a  few  of  them, 
because  they  can  be  found  only  by  search  procedures. 
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In  che  study  of  the  earth's  crust,  seismologists  face  eigenvalue 
problems  of  this  type  whpr.  they  calculate  dispersion  curves  for  Love 
waves  in  a  layered  medium  [16].  Dispersion  curves  show  the  relation 
between  the  frequency  or  period  of  a  harmonic  wave  and  its  phase 
velocity  or  ’wave  number.  Only  the  fundamental  wave  mode,  which  is  the 
mode  of  greatest  phase  velocity  or  smallest  wave  number  at  any  given 
frequency,  and  one  or  two  higher  modes  are  usually  considered.  Haskell 
[19]  has  formulated  a  matrix  transfer  method  suited  for  automatic 
computation  that  permits  finding  a  few  of  the  lower  wave  modes  in  a 
layered  elastic  medium. 


Expansion  of  Displacements  into  Eigenfunctions 


The  actual  displacements  in  the  layered  region  R  shown  in  Fig.  4 
may  be  expanded  into  eigenfunctions,  which  are  determined  from  the 
solutions  to  Eq.  (4.16)  and  which  satisfy  the  differential  equations, 

Eq.  (4.2),  and  the  boundary  conditions  at  horizontal  planes,  Eq.  (4.11). 
According  to  Eq.  (4.9),  the  expansion  is  of  the  form 


S(x,z,t) 


a  v  (z)  exp(iwt-ik  x) 
s  s  s 


S=1 


(4.17) 


in  which  vs(z),  kg  and  as  are  the  mode  shape,  the  wave  number  and  the 
mode  participation  factor  of  the  sth  mode  or  eigenfunction,  respectively. 

The  mode  participation  factors  have  to  be  chosen  so  that  dis¬ 
placement  and  stress  conditions  at  the  boundary  between  the  irregular 
and  the  layered  region  are  satisfied.  In  addition,  only  those  modes 
which  transmit  energy  in  the  positive  x-direction  and  do  not  increase 
in  amplitude  with  x  can  oe  included  in  the  expansion.  These  latter 
conditions  follow  from  energy  considerations.  Since  the  motion  is 
generated  by  external  forces  in  the  irregular  region  and  since  the 
layered  region  R  is  open  to  the  right,  the  energy  transmission  must  be 
positive  in  the  x-direction  and  the  energy  density  cannot  increase  with  x. 
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The  method  outlined  above,  which  is  based  on  continuum  theory, 
was  employed  to  solve  a  very  simple  boundary  value  problem  involving 
a  line  load  on  a  homogeneous  layer,  see  Appendix  3.  However,  major 
difficulties  arise  when  a  more  complicated  problem  is  to  be  solved  by 
this  method.  In  actual  analyses,  the  number  of  eigenfunctions  included 
in  the  expansion,  Eq.  (4.17),  will  be  small,  instead  of  infinite,  because 
the  determination  of  any  one  eigenfunction  takes  a  considerable  effort. 

In  addition,  it  is  difficult  to  ensure  that  none  of  the  more  important 
eigenfunctions  is  omitted  in  the  expansion  and  that  a  good  approximate 
solution  is  obtained. 

Therefore,  a  discrete  method  will  be  developed  for  the  analysis 
of  harmonic  motion  in  the  layered  region.  This  method  is  suited  to 
automatic  computation,  provides  accurate  results  and  is  easy  to  use  in 
connection  with  the  finite  element  analysis  of  the  irregular  region. 


4.2  Eigenvalue  Problem  for  Discretized  Layered  Region 

The  layered  region  is  discretized  in  the  vertical  direction  by 
subdividing  the  region  into  thin  layers  as  shown  in  Fig.  11  and  assuming 
that  v(z)  varies  linearly  within  each  layer.  In  order  for  this  assump¬ 
tion  to  be  reasonable  the  thickness  ,h,  of  each  layer  must  be  chosen 
small  compared  to  the  wavelength  of  shear  waves  in  the  iyer.  For  this 
reason  the  number  of  layers,  n,  in  the  discrete  system  is  usually  larger 
than  the  number  of  natural  layers;  a  typical  analysis  involves  from  10 
to  40  layers.  The  region  is  subdivided  into  layers  so  that  the  inter¬ 
faces  between  the  layers  coincide  with  the  nodal  joints  at  the  vertical 
boundary  with  the  irregular  region.  The  layered  and  the  irregular 
regions  can  then  be  connected  at  these  nodal  joints.  As  the  displace¬ 
ments  at  the  vertical  boundary  vary  linearly  between  the  nodal  joints 
of  the  finite  elements  and  between  the  layer  interfaces,  the  displace¬ 
ments  will  be  continuous  across  the  boundary. 

The  assumption  of  linear  variation  of  v(z)  within  each  layer 
implies  that  v(z)  is  defined  by  its  values  v^=v(z..),  j=l,  ...»  n,  at 
the  layer  interfaces.  Therefore,  the  displacements  of  a  generalized 
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Love  wave  may  be  represented  by 

{£}  =  a{v}  exp(iiot-ikx) 


(4.18) 


in  which  the  vector  {6}  contains  the  displacements  S  ^  ,  j=l,  ....  n,  of 
the  layer  interfaces.  The  vector  {v}  represents  the  mode  shape  and  a 


is  the  mode  participation  factor  as  before. 
Lsplacement  function  for 

u(x,z)  =  v(z)  exp(-ikx) 


The  displacement  function  for  z^<z<z^+^  is  accordingly 


(4.19) 


in  which 


v<_;  =  (^-o/hj  •  v.  +  (2-z.)/h.  •  v.+1 


(4.20) 


Omitting  the  common  factor  exp(iodt),  these  displacements  produce  the 
shear  strains 

3u 


c  =  =  -  iku 

xy  ox 


3u 


-  ^  =  (-v.  +  v  )/h.  *  exp(-ikx) 

zy  3z  j  j+1  j  K 


(4.21a) 

(4.21b) 


shear  stresses 
T 

T 

and  inertia  forces 
f 


=  G.  c 
xy  j  xy 

(4.22a) 

=  G.  z 
zy  j  zy 

(4.22b) 

2 

T  °j  J  U 


(4.23) 


If  a  section  of  the  layered  region  between  the  planes  x=0  and 
x=^  is  considered  separately,  see  Fig.  11,  the  surface  tractions 


and 


Po(2)  =  ~  Txy(°’z) 


pr(z)  =  +  T  (£,Z) 
s 


(4.24a) 

(4.24b) 


must  be  applied  at  the  boundaries  x=0  and  x=£  in  order  to  preserve 
dynamic  equilibrium.  They  are  the  only  external  forces  acting  on  this 
section. 

The  application  of  the  principle  of  virtual  work  as  expressed  by 
Eq.  (3.26)  yields 
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(4.25) 


r*  rZ  n+1  rZ  n+1 

/  /  (e*  T  +  £*  T  -  u*fT)dz  d x  =  /  u*(p  +  Pr)dz 

JQ  J7  xy  xy  zy  zy  I  Jz ^  *o  QZ 

in  which  £*  ,  e*  and  u*  are  the  complex  conjugates  of  the  virtual 
xy  zy 

strains  and  displacements,  respectively. 

According  to  Eq.  (4.18),  the  only  degrees  of  freedom  in  the  dis¬ 
cretized  region  are  the  nodal  displacements  v. ,  j=l,  ...»  n.  From  the 
virtual  displacement  state  "j",  i.e. 


v.  =  1  and  v^  =  0  for  £.  ^  j 


(4.26) 


and  from  Eqs.  (4.19  -  4.21)  it  follows  that 


(z-z  )/h  *exp(ikx) 

for  z.  <z<z. 

1-1 - 1 

U*  =  < 

( z j-z) /h j  *  exP ( ikx) 

for  z  ,<z<z . . , 

3-  ~  3+1 

(4.27a) 

0 

elsewhere 

e*  = 
xy 

ik  u* 

(4.27b) 

and 

1/tu  ^*exp(ikx) 

for  z.  <z<z. 

3-1  3 

£*  =  • 
zy 

-1/h^  *exp(ikx) 

for  z.<z<z... 

3  3+1 

(4.27c) 

0 

k 

elsewhere 

in  which 

k  is  the  complex  conjugate 

of  k. 

Substitution  of  Eqs.  (4.19  -  4.24)  and  (4.27)  into  Eq.  (4.25)  leads 
to  the  equilibrium  equations  for  the  section  considered.  The  n  virtual 
displacement  states  "j",  j=l,  2,  ...,  n,  produce  a  set  of  simultaneous 
linear  equations  in  the  n  unknowns  v. .  When  the  simple  integrations 
defined  in  Eq.  (4.25)  are  performed,  the  equilibrium  equations  can  be 
reduced  to  the  matrix  equation 

(k2[A]  +  [G]  -  w2[M]){v}  =  {0}  (4.28) 

in  which  {v}  is  the  displacement  vector  with  the  elements  v^ ,  j=l,  ...,  n. 
The  n  x  n  matrices  [A],  [G] ,  and  [M]  consist  of  the  contributions  from 
the  individual  layers  and  can  therefore  be  conveniently  assembled  from 


the  layer  submatrices  as  demonstrated  in  Fig.  13.  The  submatrices  to 
be  substituted  for  [X]j,  j=l,  ...,  n,  in  Fig.  13  are 


G. 

CGli  "  h1 
1 


1MiJ  “  Vi 


r  i  i 

3  6 

1  1 
-  A  1 


-1  1 


j— 1 j  • * .  »  n 


(4.29a) 


j~lj  .  •  •  »  n 


j “1 »  ...»  n 


(4.29b) 


(4.29c) 


for  the  matrices  [A],  [G] ,  and  [M] ,  respectively.  The  first  two 
matrices  are  obviously  related  to  the  stiffness  of  the  layers.  They  :.j.e 
real  in  the  undamped  case,  when  all  shear  moduli  G__,  j=l,  ...»  '  "re 
real.  The  mass  matrix  [M]  is  always  real  since  ,  the  mass  density  of 
the  jth  layer,  is  real. 

Eq.  (4.28)  is  independent  of  £.  Hence,  the  solutions  to  Eq. 

(4.28)  satisfy  equilibrium  in  the  layered  region  between  the  vertical 
planes  x=0  and  x=£  where  £  may  take  any  value  greater  than  zero.  In  the 
following  it  will  be  assumed  that  £  is  infinitely  large. 

As  the  circular  frequency  -jj  is  a  given  parameter,  it  is  convenient 


to  introduce  the  n  x  n  matrix 


[C]  =  [G]  -  uT[M] 
and  write  Eq.  (4.28)  in  the  form 

([A]k2  +  [CD  {v}  =  {0} 


(4.30) 


(4.31) 


which  states  an  algebraic  eigenvalue  problem  with  n  eigenvalues  k“,  s=l, 
...»  n,  and  the  corresponding  eigenvectors  (v}g. 

The  matrices  [A]  and  [Cj  are  tridiagonal  and  symmetric.  In  the 
undamped  case  they  are  real  and  matrix  [A]  is  positive  definite,  giving 
F.q.  (4.31)  the  convenient  property  that"  all  the  eigenvalues  and  eigen¬ 
vectors  are  real  [49]. 


37 


A  procedure  for  computing  the  complete  set  of  eigenvalues  and 

corresponding  vectors  in  the  damped  and  undamped  case  is  presented  in 

Appendix  1  and  a  FORTRAN  IV  subroutine  for  automatic  computation  is 

listed  in  Appendix  5.  The  procedure  uses  Newton's  method  to  find  the 

2 

roots  of  the  determinant  |[A]k  +  [Cj |  and  inverse  iteration  to  obtain 
the  eigenvectors. 

The  orthogonality  conditions  of  the  eigenvectors  can  be  derived 
from  any  two  solutions  to  Eq.  (4.31)  such  as 

([A]k*  +  [C]){v>r  =  {0}  (4.32a) 

and 

([A]kg  +  [C]){v}g  =  {0}  (4.32b) 

The  first  equation  is  premultiplied  by  {v}g,  the  second  equation  is 
transposed,  post-multiplied  by  {v}r  and  then  subtracted  from  the  first 
equation.  This  gives 

(kj-kg){v>*  [A] {v>r  =  0  (4.33) 

which  implies 

T  1  for  r=s 

{vriA]{v}  =  (4.34) 

0  for  r^s 

when  the  case  r=s  is  used  to  normalize  the  vectors. 

The  general  motion  of  the  discretized  layered  region  can  be  des¬ 
cribed  by  a  linear  combination  of  Love  modes.  Hence,  if  only  one  each 
of  the  modes  corresponding  to  the  pairs  iks>  s=l,  ...»  n,  exist  in  the 
system,  the  displacements  at  x=0  are 

n 

(u}R  =  \  a  {v}  =  [V]{a}  (4.35) 

Z_,  s  s 

s=l 

in  which  {a}  is  a  column  vector  containing  the  generally  complex 

participation  factors  ag,  s=l,  ...,  n,  and  [V]  is  an  n  x  n  matrix  which 

contains  the  mode  shapes  {v}  in  its  columns. 

s 

The  orthogonality  properties,  Eq.  (4.34),  imply  that 

[V]"1  =  [V]T  [A]  (4.36) 
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4. 3  Wave  Types 

2 

Each  eigenvalue  k  yields  two  possible  wave  numbers,  +  k  and  -  k  • 
s  s  s 

One  sign  corresponds  to  a  generalized  Love  wa\ e  which  propagates  or 

decays  to  the  left,  the  other  sign  to  a  wave  which  propagates  or  decays 

to  the  right.  Both  waves  have  the  same  mode  shape. 

In  the  undamped  case  all  eigenvalues  and  mode  shapes  are  real,  as 

the  matrices  [A]  and  [C]  are  real.  Hence,  by  Eq.  (4.18),  positive 
2 

eigenvalues,  k  >0,  correspond  to  propagating  waves  of  the  type 

{5}  =  a  {v}  exp(iwt  ±  i|k  |x)  (4.38a) 

s  s  s  s 

2 

Negative  eigenvalues,  ks<0,  yield  purely  imaginary  wave  numbers 
and  motions  of  the  type 

{6}  =  a  {v}  exp(iajt  ±  jk  lx)  (4.38b) 

s  s  s  s 

which  decay  or  increase  exponentially  in  amplitude  with  x  and  do  not 
propagate. 

in  the  damped  case,  when  the  matrices  [A]  and  [C]  are  complex,  all 

eigenvalues  and  therefore  all  wave  numbers  are  complex,  i.e. 
s 

k  =k~  +  ik„.  The  motion  is  a  wave 
si  2 

{5}s  =  as^vJs  exp(iwt-ik®x+k*x)  (4.38c) 


g 

which  propagates  in  the  x-direction  with  the  phase  velocity  cs=o)/k^  and 

g 

a  decaying  or  increasing  amplitude,  depending  on  the  sign  of  k^. 

In  the  undamped  case  the  possibility  exists  that  ks=0  and  that 
the  motion  degenerates  to 

{6}  =  a  {v}  exp(iwt)  (4.38d) 

s  s  s 

This  motion  is  independent  of  x  and  consists  of  a  shear  wave  travelling 
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vertically  up  and  down  through  the  layers,  being  reflected  at  the  rigid 
base  and  the  free  surface.  This  wave  type  can  only  occur  at  certain 
frequencies  which  are  the  natural  frequencies  for  one-dimensional 
vertical  waves  in  the  layered  system.  The  natural  frequencies  can  be 
found  from  the  secular  equation 

| [G]  -  w2[M] |  «  0  (4.39) 

which  follows  from  Eqs.  (4.30)  and  (4.31). 


4.4  Dynamic  Stiffness  Matrix  for  Layered  Region 


In  order  to  develop  the  dynamic  stiffness  matrix  for  the  layered 
region  R  in  Fig.  11  it  is  first  assumed  that  only  those  modes  which 
decay  or  propagate  energy  in  the  positive  x-direction  exist  in  R.  This 
requires  choosing  those  n  wave  numbers  out  of  the  2n  numbers 
±  ks=±(k^+i  k^K  s=l,  ...»  n,  that  have  a  negative  imaginary  part  or,  if 
the  imaginary  part  is  zero,  have  a  positive  real  part,  i.e.  for  s=l,  ..., 
n 


kl"i 1 k2 1 

if 

k2*° 

.  +  |kj| 

if 

(4.40) 


The  displacements  in  the  layered  region  are  then  expressed  by  super¬ 
position  of  the  corresponding  modes 

n 


{6} 


2 

S=1 


{v)sexp(iut-iksx)  ag 


(4.41) 


The  only  stress  on  the  plane  x=0,  see  Fig.  15,  is  the  shear  stress 


T  =  G  (x=0)  =  -  i  G  \  k  v  (z)  a  (4.42) 

xy  ox  Z*  s  s  s 

s=l 


in  which  G  is  the  shiar  modulus.  The  time  factor  exp(iu>t)  is  understood. 
The  assumption  that  v(z)  varies  linearly  within  a  layer  implies  that 
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also  varies  linearly  with  z  as  shown  m  Fig.  15a.  The  nodal  forces 
which  are  in  equilibrium  with  the  stresses  at  rne  vertical  face  of  the 
jth  layer  are 


P!  =  7  G.h.  Y  k  (2  vS  +  v®. .)  a 
3  6  3  3  sL,  s  3  3+1  s 


(4.43a) 


P"  =  7  G.h.  y  k  (v®  +  2  v®  )  ct 
3+1  6  3  :  <6  S  3  3+r  s 


(4.43b) 


in  which  v®  is  the  j ch  element  of  the  eigenvector  {v}s. 

These  forces,  which  are  also  shown  in  Fig.  15a,  act  on  the  region 
xN0.  The  total  forces  acting  at  the  boundary  nodal  joints  of  the 
layered  region  at  x=0  are 


P  =  P’  +  PV 
3  J  i 


j=l,  ...»  n 


with  P-pO.  These  forces  can  be  expressed  in  matrix  notation  by 
(p}R  =  i [A]  [V]  [K]  {ct} 


(4.44) 


(4.45) 


in  which  (P }  contains  the  forces  P..  3=1,  n,  and  [K]  is  a  diagonal 

matrix  with  the  elements  ks,  s=l,  ....  n,  which  3re  chosen  according  to 
Eq.  (4. 40).  Matrix  [Vj  contains  the  eigenvectors  {v;s,  s=l,  ....  n, 
columnwise  end  [A]  is  the  matrix  defined  above  by  Eq.  (4.29a). 

The  boundary  forces  {p}  are  expressed  in  Eq.  (4.45)  in  terms  of 
the  mode  participation  factors  ct^.  However,  in  the  finite  element 
analysis  of  the  irregular  region  the  nodal  displacements  are  introduced 
as  the  unknowns.  The  coordinate  transformation  according  to  Eq.  f4.27) 


yields 


in  which 


{P}R  =  [R](u}R 


m  -  i[A]  [V]  [K  j  IV]T [A] 


(4.46) 


(4.47) 


is  the  dynamic  stiffness  matrix  of  the.  semi-infinite  layered  region. 

The  matrix  IR]  relates  the  r.odal  forces  acting  cn  the  layered  region  to 
the  simultaneous  nodal  disolncements  at  the  boundary  x=0.  The  attribute 
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"dynamic"  indicates  that  the  effect  of  the  inertia  forces  is  included  in 
[R].  Matrix  [R]  is  symmetric,  as  it  should  be  according  to  the  theorem 
of  reciprocity  [14].  Its  elements  are  generally  complex.  The  real  part 
of  [R]  exactly  represents  the  elastic  resistance  of  the  discretized 
semi- infinite  layered  region  to  the  displacement  of  the  nodal  joints. 
Whereas  the  imaginary  part  of  [R]  represents  the  damping  effect  of  the 
layered  region  on  the  displacements  of  the  nodal  joints.  This  damping 
is  caused  by  dissipation  of  energy  to  infinity  or  by  viscosity  of  the 
material. 

The  dynamic  stiffness  matrix  changes  with  the  frequency  of  the 
harmonic  motions  because  the  displacement  functions  chosen  are  the 
frequency  dependent  eigenfunctions  of  the  layered  system.  Eq.  (4.46)  is 
also  valid  for  the  static  case,  co=0. 

The  dynamic  stiffness  matrix  becomes  singular  if  any  one  of  the 
wave  numbers  kg  of  the  diagonal  matrix  [K]  vanishes.  If  this  is  the 
case  no  forces  are  required  to  excite  the  mode  corresponding  to  kg-0. 
However,  this  can  happen  only  in  the  undamped  case  and  then  only  at  the 
natural  frequencies  determined  by  Eq.  (4.39). 

The  analysis  of  a  left  layered  region  L,  as  shown  in  Fig.  4,  is 
analogous  to  that  of  a  right  layered  region.  The  only  difference  between 
a  left  and  a  similar  right  layered  region  is  due  to  their  positions  with 
respect  to  the  global  coordinate  system  introduced  in  Fig.  4.  Thus  the 
dynamic  stiffness  matrix,  [L],  of  a  left  layered  region  is  also  given 
by  the  right  hand  side  of  Eq.  (4.47). 


4.5  Energy  Transmission 

Consider  a  set  of  harmonic  external  forces 

(f)  =  {P}  exp(iut)  =  ({P^}  +  i{P2))  exp(iwt)  (4.48) 

which  act  on  the  nodal  points  of  a  discrete  system  and  produce  the  nodal 
displacements 

{6}  =  {u}  exp(iwt)  =  ({u,}  +  i{u2))  exp(icot) 
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(4.49) 


r?  ■£,-wt*t^ 


The  time  average  of  the  rate  of  work  done  by  the  forces  on  the  displ 


ace- 


raeuf .'  is 


T 

E  =  jp  f  Re  {6}T  Re  {F}  dt 
0 


(4.50) 


in  which  T=2ir/u)  is  the  period  and  {6}  is  the  time  derivative  of  {6}.  Sub¬ 
stitution  of 


and 


Re  {6}  =  -  ai( { u. }  sinwt  +  {u„}  cosuit) 
1  £, 


Re  {F}  =  {Pj,}  cosuit  -  {P^}  sinuit 
into  Eq.  (4.50)  and  integration  over  one  period  yields 

E  =  §  ({u1}T{P2}  -  {u2)T{ P^) 
which  is  identical  to 

E  -  |  Im  ({u}*{P}) 

The  *  indicates  the  conjugate  transposed  and  Im  denotes  the 


(4.51) 

(4.52) 


(4.53) 


(4.54) 

imaginary 


part 


The  energy  transmission  from  the  irregular  region  I  to  the  layered 
region  R  at  the  boundary  x=0  is  now  considered,  see  Fig.  13.  Suhstitution 
cf  Eqs.  (4.35)  and  (4.45)  into  Eq.  (4.54)  gives  the  time  average 

E  =  |  Im  (i{a)*[VJ*[A] [V] [K] {a}) 


(4.55) 


This  expression  can  be  simplified  in  the  undamped  case,  when  [A]  and 
therefore  [V]  are  real,  to 


E  =  |  Im  (i{a}*[K]{a})  (4.56) 

since  [V] 4  f A] [V]  =  [I]  according  to  Eq.  (4.34).  As  [K]  is  a  diagonal 


ma 


trix,  Eq.  (4.56)  may  also  be  written  in  the  form 


««  u 

1  ’  ‘i  1  <WKI2  - 1  E= 

s=l 


(4.57) 


s=l 
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5.  AXI SYMMETRIC  LOVE  WAVE  MOTION 


5. 1  Motion  in  Layered  Region 

The  analysis  of  axisymmetric  Love  wave  motion  in  a  layered  region 
as  shown  in  Fig-  5  is  similar  to  the  analysis  of  plane  Love  wave  motion 
presented  in  Chapter  4. 

Any  point  (r,z)  in  the  region  undergoes  displacements  in  the 
angular  9-direction  only.  The  displacements  of  harmonic  motion  at  the 
frequency  o>, 


6  =  6A(r,z,t)  =  u(r,z)  exp(iwt) 


(5.1) 


are  governed,  according  to  Eq.  (3.16c),  by  the  homogeneous  differential 
equations  for  the  n  layers 


o  /lit  .  1  3u  u  5  u  \  2  _  n 

j  U  2  r  3r  2  .2' 


j=l,  .. . ,  n  (5.2) 


The  assumption  that 


u(r,z)  =  v(z)  •  g(r) 


leads  to 


(5.3) 


v  G.  g  \  -2  r  dr  2 / 


j=l>  • . . ,  n  (5.4) 


which  can  only  be  valid  for  all  values  of  r  and  z  in  the  respective 

layers  if  both  sides  ol  the  equation  are  equal  to  the  same  constant. 

2 

Setting  the  constant  equal  to  k  ,  the  following  ordinary  differential 
equations 


2 

+  (p.O)2/G.-k2)v  =  0 
dz2  3  3 


j=l,  . . . ,  n  (5.5) 


and 


.2 


^-f  +  -  4s  -  (-4  -  k2)g  =  0 

,2  r  dr  2 
dr  r 


(5.6) 


are  obtained.  Kq.  (5.6)  is  a  Bessel  equation  [2]  which  has  a  solution 

(5.7) 


g(r)  =  H^2) (kr) 


,(2) 


in  which  H^-''  is  the  Hankel  function  of  the  first  order  and  second  kind 


[2].  Eq.  (5.7)  can  easily  be  verified  by  substitution  into  Eq.  (5.6) 
and  observing  the  following  derivatives 


£  »‘2)o-> 


-  k  (kr) 


(5.8a) 


~  H  j2)  (kr)  =  -  -  H(2)(kr)  +  k  H(2)(kr) 
dr  1  rl  o 


(5.8b) 


-4  ^  i 2>  (kr)  -  (—  -  k2)  h|2)  (kr)  -  £  H(2)(kr) 
,21  2  1  r  o 

dr  r 


(5.8c) 


m  (2) 

The  Hankel  functions  (kr)  and  (kr)  tend  to  zero  in  the 


sector  -tt<  arg  (kr)  <0  as  ]kr|  -*».  They  are  related  to  the  Bessel 


functions  of  the  first  and  second  kind,  J  and  Y  ,  respectively,  by 

m  in 


H(2)  =  J  -  i  Y 
m  m  m 


(5.9) 


in  which  the  subscript  ra  indicates  the  order  of  the  functions  and 


i  =  /-T.  However,  Eq.  (5.9)  should  not  be  used  to  compute  values  of 

,(2) 


H'“/(kr)  in  the  sector  -u<  arg  (kr)  <0  from  values  of  J  and  Y  because 
m  mm 


severe  cancellations  may  occur  in  the  subtraction  J  -  i  Y  and  may 

J  mm 


result  in  the  loss  of  all  significant  digits.  A  method  for  computing 


(2)  (2)  ... 

the  functions  Hq  and  by  ascending  and  asymptotic  series  is 


presented  in  Appendix  4  and  a  Fortran  IV  subroutine  for  automatic  com¬ 
putation  of  these  functions  is  listed  in  Appendix  5. 

A  solution  to  Eq.  (5.2)  may  be  expressed  in  the  form  of  an  axi- 
symmetric  wave 

/ON  /9\ 

(5.10) 


6  =  v(z)  (kr)/H^  (krQ)  •  exp(icot) 


A6 


t . ,  i -iinirrniiffiaiai^  . ^ 


feaF??^4--4 


which  is  similar  to  a  plane  wave  at  some  distance  from  the  origin.  In 


this  expression  rQ  defines  the  distance  of  the  boundary  between  the 


irregular  and  layered  regions  from  the  origin.  If  |kr|  is  large  the 
asymptotic  approximation  to  the  Hankel  function  [2] 

,(2) 


H,  ;  (kr)  =  /2/  (krTT)  exp(-ikr+i3Ti/4) 


(5.11) 


may  be  used  to  obtain 


6  =  v(z)  exp(iwt-ikr)  /2/ (trkr)  exp(i3iT/4)/H 


(2) 


(kr  ) 

o 


(5.12) 


which  is  similar  to  Eq.  (4.9)  except  for  the  decay  factor  /2/ (irkr)  and 
a  complex  number  indicating  a  phase  shift.  The  similarity  to  a  plane 
wave  shows  that  Eq.  (5.10)  describes  a  wave  travelling  away  from  the 
origin  if  the  real  part  of  the  wave  number  k  is  positive. 


The  boundary  conditions  at  the  horizontal  planes  z. ,  j=l,2,  ...» 
n+1,  are,  as  in  the  plane  case,  the  continuity  of  the  displacements  and 


the  stresses  at  the  layer  interfaces,  zero  shear  stresses  at  the  free 
surface  and  zero  displacements  at  the  fixed  base.  As  the  stress- 
displacement  relation  is 

Tz0  =  G  Iz  =  G  dz'  HP^(kr)/Hi2^kr0)  *  exp(iwt)  (5.13) 


the  boundary  conditions  take  the  form 


and 


V(z')  =  V(zp 

j_2 j  •  ♦  •  j  n 

(5.14a) 

Gdv(zP  =  Vi  to v( 

j=2 9  i  n 

(5.14b) 

v(Vi}  =0 

1 

(5.14c) 

iz  v(zP  - 0 

(5. 14d) 

The  function  v(z)  is  determined  by  the  n  ordinary  differential 
equations,  Eq.  (5.5),  and  the  boundary  conditions,  Eq.  (5.14),  which  are 
identical  to  the  corresponding  differential  equations  and  boundary  con¬ 
ditions  for  the  plane  case,  Eqs.  (4.5)  and  (4.11).  Therefore,  the  re¬ 
sulting  eigenvalue  problem  and  its  solutions  are  the  same  for  the  plane 
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and  axisymmetric  cases.  This  holds  not  only  for  the  continuum  theory 
but  also  for  the  discrete  theory  because  in  the  axisymmetric  case  the 
layered  region  is  discretized  in  the  same  manner  as  in  the  plane  case. 

It  is  again  assumed  that  v(z)  varies  linearly  within  each  layer 
according  to  Eq.  (4.20).  Hence,  the  displacements  in  the  layered  region 
may  be  defined  by  the  displacements  of  the  layer  interfaces  which  are 
collected  in  the  vector  {5},  i.e. 


T 

{5}  =  <5.6.  ...  6  > 

12  n 


(5.15) 


The  sth  axisymmetric  mode  is  obtained  from  the  sth  solution  to 
Eqs.  (4.31)  and,  by  Eq.  (5.10),  takes  the  form 


{5}  =  {v}  Hp^(k  r)/Hf“^(k,r  )  •  exp(itot) 

s  si  s  i  lo 


(5.16) 


Thus,  the  displacements  in  the  layered  region  may  be  expressed  by  the 
linear  combination 


{c}  =  *)>  a  {v}  nf (k  r)/Hp\k  r  )  •  exp(iwt)  (5.17) 

/_  ^  S  S  -L  S  «L  SO 


The  rule  for  the  selection  of  the  n  wave  numbers  k  from  the  2n  values 

s 

ik  ,  s=l,  ...,  n,  is  the  same  as  in  the  plane  case,  i.e.  Eq.  (4.40). 

S  (2^ 

This  follows  from  the  properties  of  '  and  its  similarity  to  the 

exponential  function  for  large  arguments. 


5.2  Dynamic  Mi. ffness  Matrix  for  Layered  Region 


The  nodal  displacement  at  the  cylindrical  boundary  r=ro  between  the 
layered  and  irregular  regions  are,  omitting  the  time  factor  exp(iiot). 


(u}R  =  \  (v) ji  =  fV]{a} 
s  s 

s=I 


(5.18) 


in  which  the  vector  { u}  has  the  elements  u^ ,  j=l,  ...,  n,  the  n  x  n 
matrix  {V]  contains  the  eigenvectors  {v}g  columnwise,  and  the  vector 
{a}  contains  the  mode  participation  factors  ag. 
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(5.19) 


The  inversion  of  Eq.  (5.18)  gives,  according  to  Eq.  (4.36), 

{a}  =  IV]T[A]{u}R 

which  serves  to  obtain  the  mode  participation  factors  from  a  known 
vector  (u‘  . 

The  only  stress  at  the  boundary  r=r0  is  the  shear  stress  Trg.  By 
using  the  stress-displacement  relation 


■i-  =  G  (~  -  -) 

ro  ur  r 


(5.20) 


and  Eq.  (5.17),  the  shear  stress  in  the  jth  layer  at  r=rQ  is  obtained. 


t  a  =  G  y  v  (z)  (k  H(2)(k  r  )/H^(k  r  )  - 
ro  s  so  so  l  so 


— }  a 
r  s 
o 


(5.21) 


when  the  time  factor  exp(iu)t)  is  omitted.  The  assumption  that  v(z)  varies 
linearly  within  a  layer  implies  that  Trg  also  varies  linearly  with  z  as 
shown  in  Fig.  15a.  The  nodal  forces  which  are  equivalent  to  the  stresses 
acting  on  the  vertical  boundary  of  a  one  radian  segment  of  the  jth  layer 


G.h./6  -y  (2vS  +  v®  ){-k  r  •  H(2)(k  r  )/H7(2)(k  r  )  +  2}  a 
J  3  3  3+1  soo  sol  so  s 

s=l 

(5.22) 


P"  =  G.h./6  (vS  +  2v®. ){-k  r  •  H(2)(k  r  )/H.(2)(k  r  )  +  2}  a 
J+1  j  j  <L>  j  j+1  so  o'so'l'so  s 


(5.23) 


in  which  v.  is  the  jth  element  of  the  eigenvector  {v}  .  These  forces, 
J  s 

which  are  also  shown  in  Fig.  15a,  act  on  the  region  t>rQ. 

The  total  forces  acting  on  a  one  radian  segment  of  the  layered 
region  are 


P  -  P’  +  P" 

J  j  j 


j  1 ,  •»•»  ti 


(5.24) 
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,•  J.V™ 


with  P^'=0.  In  matrix  notation  these  forces  may  be  expressed  by 

{P)R  =  [A] [V] [H ] {a}  (5.25) 

The  vector  { P }  contains  the  forces  P^.,  j=l,  •*•»  n;  the  n  x  n  matrix 
[V]  contains  the  eigenvectors  {v}s,  s=l,  •••»  n,  columnwise;  [A]  is  the 
matrix  defined  by  Eq.  (4.29)  and  [H]  is  the  diagonal  matrix 

[H]  =  Diag  {-k  r  H(2)(k  r  )/H.(2)(k  r  )  +  2}  (5.26) 

soo  so  1  so 

Substitution  of  Eq.  (5.19)  into  (5.22)  gives 

{P)R  =  [R]{u}R  (5.27) 


in  which 

[R]  =  [A][V][H][V]T[A]  (5.28) 

The  matrix  [R]  is  the  dynamic  stiffness  matrix  for  a  one  radian 
segment  of  the  axisymmetric  layered  region  which  extends  in  tne  radial 
direction  from  rQ  to  infinity.  It  relates  the  nodal  forces  to  the 
simultaneous  nodal  displacements  at  the  boundary  r=ro  and  is  valid  for 
axisymmetric  harmonic  motion  including  the  static  case.  Matrix  [R] 
changes  with  frequency  and  depends  on  rQ.  It  cannot  become  singular, 
because  none  of  the  diagonal  elements  of  [H]  can  vanish  at  any  given 
frequency  and  matrices  [A]  and  [V]  are  always  nori-singular. 


y.JW,.—.  jtoJgsr^j^iife^r-K-^' Vis' ^E^JKSJ-^^sSjg^ 


6.  PLANE  RAYLEIGH  WAVE  MOTION 


6 . 1  Eigenvalue  Problem 

Free  harmonic  Rayleigh  wave  motion  under  p] nne  strain  conditions  in 
a  semi-infinite  layered  region  as  shown  in  Fig.  5  is  considered.  Any 
point  in  the  region  undergoes  displacements  in  the  x  and  z-directions , 


6  =  q  (x, z)  exp(iwt) 

X  X 

<$z  -  qz(x,z)  exp(ioit) 


(6.1a) 

(6.1b) 


The  homogeneous  differential  equations  of  motion  for  the  spatial 
parts  of  the  displacements  are  according  to  Eq.  (3.15),  for  the  jth  layer. 


/  32q  32q  \  /  32q  32q  \  „ 

ci  \jf  *  777  +  (7 +  V  W +  rai)  +  PJ  “  "*  ■ 0  (6-2a) 


/32q  32q  \  /  32q  32q  \ 

G  (— ~  +  +  (X.  +  G.)  (  +  — T-)  +  p.  0)2q  =  0  (6.2b) 

J  \3x2  3z2  /  J  J  \8xoz  3z2  /  J  z 


Separation  of  the  variables  x  and  z  with  the  assumptions  that 


q  (x,z)  =  u(z)  g (x) 


(6.3a) 


q  (x,z)  =  w(z)  g(x) 


(6.3b) 


leads  to  the  coupled  ordinary  differential  equations 

2 

(k2U.  +  2G.)  -  aj2p.}  u-G.  ~  +  ik  (X.  +  G.)  —  =  0 

J  J  J  J  dz2  J  J  dz 


(6.4a) 


,2 

aw  ,  .. 


{k  G.  -  up.}  w  -  (X.  +  2G.)  2JS.  +  ik  (X.  +  G.)  ~  =  0 

J  J  J  J  .2  J  3  dz 


(6.4b) 
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(Xi  +  2Gi'  dl  w(*i)  "  ik  xi  u(zi}  “  0  (6.9e) 

and 

-ik  w(zp  +  ^  u(z'p  =  0  (6. 9f) 

in  which  zj  and  z'j  respectively  refer  to  planes  immediately  above  and 
immediately  below  the  interface  z.. 

Eqs.  (6-4)  and  (6.9)  lead  to  an  eigenvalue  problem  similar  to  that 
for  Love  waves.  In  continuum  theory, this  problem  consists  of  2n  simul¬ 
taneous  homogeneous  equations  with  coefficients  which  contain  the  eigen- 
2 

value  k  in  the  argument  of  transcendental  functions.  A  derivation  of 
the  eigenvalue  problem  can  be  found  in  Ref.  [22].  Finding  solutions  to 
this  problem  is  even  more  difficult  than  finding  solutions  to  the  corres¬ 
ponding  eigenvalue  problem  for  Love  waves,  Eq.  (4.16),  because  in  the 
Rayleigh  wave  case  th"  number  of  equations  for  n  layers  is  2n  c.r.'1  co- 
eff  cients  are  more  complicated  than  those  in  the  Love  wave  case.  For 
the  same  reasons  given  in  Section  4.1,  a  discrete  theory  is  preferred  to 
the  continuum  theory  and  will  now  be  developed. 


Discretization 


As  in  the  Love  wave  case,  the  layered  region  is  treated  as  a 
continuum  in  the  horizontal  direction  but  discretized  in  the  vertical 
direction  by  assuming  that  u(,z)  and  w(z)  in  Eq.  (6.7)  vary  linearly 
within  each  la''er.  In  order  for  this  assumption  to  be  reasonable  the 
layers  have  to  be  chosen  thin  compared  to  the  length  of  S-waves  in  the 
layer. 

The  displacements  in  the  layered  region  may  be  defined  by 


{6} 


T 


<6.6. 
xl  zl 


6  6  > 

xn  zn 


(6.10) 


in  which  6  .  and  6  j=l,  ....  n,  are  the  horizontal  and  vertical  inter- 
xj  zj 

face  displacements,  respectively.  IT  the  elements  of  the  corresponding 
mode  shape  {v}  are 


V2j-1  '  U(V 


an  a 


-j 


W(z.) 


.1=1, 


(6.11) 
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the  discrete  representation  of  a  plane  generalized  Rayleigh  wave  is 


{6}  =  (v)  exp(io)t-ikx) 


(6.12) 


The  displacement  functions  in  terms  of  interface  displacements  are 


qx(x,z)  -  u(z)  exp(-ikx) 

q  (x,z)  =  w(z)  exp(-ikx) 
z 


in  which,  for  z.<z<z.11, 
J-  ~  J+l 


(6.13a) 

(6.13b) 


u(z)  =  (Zj+]  -  z)/hj  *  v2j_i  +  (z-z/)/h.  ’  v2j+l  (6.14a) 


w(z)  =  (z.,.,  -  z)/h.  •  v_.  +  (z-z.)/h. 


j  2j 


T'  j  2j+2 


(6.14b) 


Omitting  the  common  factor  exp(i^t),  these  displacements  produce  the 


strains 


£  =  —  =  “ikq 

X  dX  X 


(6.15a) 


£z  3z 

3q,.  3<l 


=  (-v2  +  v2j4-2)  /^j  ’  exp(-ikx) 


\z  =  ST  +  W  ’  (-V2j-l  +  v2j+I)/h.-exP(-ilai)  -ikq, 


(6.15b) 


(6.15c) 


stresses 


c  =  (A.  +  2G. )  e  +  A.  e 

X  j  j  x  j  z 

a  =  (A.  +  2G.)  £  +  A.  £ 
z  j  j  z  j  x 


(6.16a) 

(6.16b) 


T  =  G.  y 

J  xz 


(6.16c) 


and  inertia  forces 

CI  2 

f  =0)  p.  q 
x  3  ^x 

CI  2 

f  =  u  p.  q 

Z  J  z 


(6.17a) 


(6.17b) 


If  a  section  of  the  layered  region  between  the  planes  x=0  and  z=£ 
is  considered  separately,  see  Fig.  12,  one  has  to  apply  as  surface 
traction  -0  (0,z)  and  -T(0,z)  at  the  left  boundary,  x=0,  and  0  (£,z)  and 

X  X 

t(£>z)  at  the  right  boundary,  x=£.  The  surface  tractions  are  the  only 


external  forces  acting  on  this  section  of  t:ie  layered  region. 

The  application  of  the  principle  of  virtual  work  as  expressed  by 
Eq .  (3.20)  yie K  ? 


A 

C*0  +  Y*  T 

z  z  xz 


dzdx 


'n+1 


-  i 


z-z.. 


-  q*(0,z) 


(0,z> 


q,;(0,z)  T(0,z) 


+  q*(C >  z)  Ox(C,  z)  +  q*(C,  z)  T(C,  z) J  dz  (6.18) 

in  which  £*,  £*,  and  Y*  are  the  complex  conjugates  of  the  virtual 

X  Z  X  z 

strains  and  q*  and  q*  are  the  complex  conjugates  of  the  virtual  displace¬ 
ment. 


The  degrees  of  freedom  in  the  discretized  region  are  v..  j=l,  ..., 
2n,  according  to  Eq.  (6.11).  The  virtual  displacement  state  "1"  defined 
by 


A  ,  A 

V.  =  1  ;  V  = 

%  m 

0  for  m^2. 

(6.19) 

yields, 

with  Eqs .  (6.13  -  6.15), 

for  odd  values  of  l. 

i.e.  1-2 j-l 

(z-z.^/h  •  exp(ikx) 

for  z.  <z<z. 

J-l-  -  J 

■■ 

(z^+^-z)/h^  *  exp(ikx) 

for  z.<z<z . , , 

J-  -  J+l 

(6.20a) 

: 

0 

elsewhere 

Y* 

xz 


i  ;  £*  -  ik  q*  ; 

e*  =  o 

X  X 

z 

l/h._j  *  exp(ikx) 

for  z.  ,<z<z. 
J-l  J 

-i/hj  •  exp(ikx) 

for  z.<z<z.., 
J  J+l 

0 

e Isewhe re 

(6.20b) 


(6.20c) 


and  for  even  values  of  SL,  i.e.  2=2  j 


(z-z.  .)/h.  ,  *  exp(ikx)  for  z.  <z<z. 

J-l  j-1  j-l~  -  j 

J3> 
N  * 

1! 

(z  -z)/h.  *  exp(ikx)  for  z.<z<z... 

J+l  J  .1-  -  J+l 

(6.20d) 

1 

L  o 

elsewhere 

0  ;  e*  =  0  ;  y* 

x  1  xz 

=  +ik  q* 

(6.20e) 

1/hj  ^  *  exp(ikx) 

for  z.  <z<z. 

J-1  J 

e*  =  i 
2 

-1/hj  *  exp(ikx) 

for  z.<z<z.,. 

J  J+l 

(6.20f) 

0 

elsewhere 

Substitution  of  Eqs. 

(6.13  •  6.17)  and  Eqs. 

(6.20)  into  Eq.  (6.18) 

leads  Co  one  equilibrium  equation  for  each  of  the  2n  virtual  displacement 
states  "2".  The  2n  equilibrium  equations  are  simultaneous  complex  linear 
equations  in  the  2n  complex  unknowns  v^,  2=1,  . ..,  2u.  When  the  simple 
integrations  defined  in  Eq.  (6.18)  are  performed,  the  set  of  equations 
may  be  written  in  matrix  notation  as 

(k2[A]  -  ik[D]  +  ik[I)]T  +  [G]  -  u)2[M]){v}  =  (C)  (6.21) 

The  vector  {v}  contains  the  complex  displacements  v.,  j=l,  ...,  2n.  The 

J 

2n  x  2n  matrices  [A] ,  [D] ,  [G]  and  [M]  consist  of  the  contributions  from 
the  individual  layers  and  can  therefore  be  assembled  by  addition  of  layer 
submatrices  as  indicated  in  Fig.  14.  The  submatrices  to  be  substituted 
for  [X]  in  Fig.  14  are 


2(2CyKV)  0  (2G  +X  J 


6  hj 


0  2G.  0  G. 

J  J 

(2Gj+Aj)  0  2(2Gj+Aj)  0 


0 


G.  0  2G .  I 
j  JJ 


j=l,  ...»  n  (6.22a) 


56 


.i . ? 


Kl 


'i 

N 


I01^ 


[M] .  =p.h . 

J  J  J  1/6 


0  -G. 


(2G.-5-X.)  0  -(2G.+X.) 

J  3  J  J 


(2G.+X.)  0  (2G.+X.) 

3  3  3  2 


0  1/6 
1/3  0 

0  1/3 

1/6  0 


j=l,  ....  n  (6.22c) 


j=l,  ...»  r.  (6.22d) 


for  the  matrices  [A],  [D] ,  [G]  and  [M] ,  respectively.  The  first  three 
matrices  are  obviously  related  to  the  stiffness  «>f  the  layers  and  are 
real  in  the  undamped  case  when  the  shear  moduli  G^  and  the  Lame's  con¬ 
stants  X ^ .  j=i,  ...»  n,  are  real.  Matrix  [M]  is  a  consistent  mass 
matrix  [*]. 

Eq.  (6.21)  is  independent  of  £.  Hence,  the  solutions  to  Eq.  (6.21) 
satisfy  equilibrium  in  the  layered  region  between  the  vertical  planes 
x=0  and  x=£  where  £  may  take  any  value  greater  than  zero.  In  the  follow¬ 
ing  it  will  be  assumed  that  £  is  infinitely  large. 

The  circular  frequency  u  is  a  given  parameter.  When  the  2n  x  2n 


matrices 


[C]  =  [G]  -  w  [M] 


(6.24) 


[B]  =  (D}1  -  ID} 

are  introduced,  Eq.  (6.23)  can  be  written 

([A]k2  +  i [B]  k  +  [C]){v}  =  0 


(6.25) 


(6.26) 


The  matrices  [A]  and  [C]  are  symmetric  and  the  matrix  [B]  is  skew 
symmetric.  [A],  [B]  and  [C]  are  real  in  the  undamped  case  and  complex 
in  the  damped  case.  Eq.  (6.26)  differs  from  a  similar  matrix  equation 
obtained  by  Lysmer  [26]  in  that  it  includes  a  consistent  mass  matrix, 
while  Lysmer  employed  a  lumped  mass  matrix.  Eq.  (6.26)  constitutes  a 
sec  of  2n  linear  homogeneous  equations  which  have  non-trivial  solutions 
{v}  only  if  the  determinant  of  the  coefficient  matrix  vanishes.  Hence, 
for  any  given  frequency  u>,  the  secular  equation 


[A]k“  +  i [B]  k  +  [C] |  =  0 


(6.27) 


defines  the  possible  wave  numbers  for  generalized  Rayleigh  waves  in  the 
layered  region. 

Because  the  determinant,  Eq.  (6.27),  is  a  polynomial  in  k  of  order 

An,  Eq.  (6.26)  states  an  algebraic  eigenvalue  problem  which  has  4n 

generally  complex  roots  k  ,  s=l,  2,  ...,  4n.  The  corresponding  solution 

s 

vectors  {v)s>  s=l,  ...,  4n  are  called  eigenvectors  or  mode  shapes.  A 
numerical  method  for  finding  all  of  the  eigenvalues  and  eigenvectors  is 
presented  in  Appendix  2  and  a  FORTRAN  IV  subroutine  for  this  task  is 
listed  in  Appendix  6-. 

I t  is  convenient  for  the  following  wave  interpretation  to  rewrite 
Eq.  (6.26)  in  the  partitioned  form 


9  I 

k  A  +  C 
_  x  x  I 

T  I 

-ik  B  , 


A  +  C 
z  z 


(6.28) 


where  the  elements  of  {v}  are  the  horizontal  displacements  v.,  j=l,  3, 

**  J 

...,  2n-l,  and  the  elements  of  {v^}  are  the  vertical  displacements 
Vj ,  j=2,  4,  ...,  2n.  The  n  x  n  submatrices  [A^],  [  A_  ] ,  [C^] ,  IC^]  and 
(B  ]  are  obtained  from  the  matrices  [A],  [C]  and  [Bj,  respectively, 
by  permuting  the  rows  and  columns  in  the  same  way  as  the  vector  elements 


are  permuted. 


6.2  Wave  Types 

Each  eigenvalue  k  and  its  corresponding  eigenvector  {v}  define 
s  s 

a  Rayleigh  wave  mode  which  can  exist  in  the  layered  region  and  has  the 
displacements 

{6}  =  a  {v}  *  exo(iu)t-ik  x)  (6.29) 

s  s  s  s 

The  nature  of  this  motion  depends  on  the  wave  number  kg.  As  in  the  Love 
wave  case  four  wave  types  can  be  distinguished. 

a)  If  k  is  complex,  k  =  k^+i  k2»  k^O,  k9?;0,  the  motion  is  of  the 

f  orm 

{5}  =  a{v}  •  exp(k2x)  •  exp(io)t-ikj>:)  (6.30a) 

which  propagates  in  the  x-direction  with  the  phase  velocity  c=t^/k^  and 
decays  or  increases  in  amplitude  depending  on  the  sign  of  k^-  The  mode 
shape  {v}  is  complex  and  generally  neither  the  horizontal  nor  the  vertical 
displacements  are  in  phase  on  vertical  planes. 

b)  If  k  is  real,  k  =  k^,  k2=0,  the  motion  is  a  wave  which  propagates 
in  the  x-direction  with  a  constant  amplitude  and  the  phase  velocity 
c=w/k^  and  is  of  the  form 

{6}  =  a{v}  exp( iut-ik^x)  (6.30b) 

Inspection  of  Eq.  (6. 28)  shows  that  for  a  real  k  the  vector  {v^.}  is 
purely  real  while  the  vector  {v^}  is  purely  imaginary.  This  means  that 
ail  horizontal  displacements  are  in  phase  at  constant  x  and  90°  out  of 
phase  with  the  vertical  displacements.  The  particle  motion  describes  an 
ellipse  with  its  major  axis  either  parallel  or  normal  to  the  surface. 

c)  If  k  is  purely  imaginary,  k=ik2,  k^O,  the  motion  is 

{5}  =  a{v)  exp(iut)  (6.30c) 

which  varies  exponentially  in  the  x-direction  and  does  not  propagate. 

As  the  coefficient  matrix  of  Eq.  (6.23)  is  real  for  a  purely  imaginary  k. 


aasalaai Ss  <  ■  • 
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and  as  the  phase  does  not  change  in  x,  the  horizontal  and  vertical  dis¬ 
placements  of  all  points  in  the  layered  region  have  the  same  phase.  In 
this  case  particles  oscillate  along  straight  lines, 
d)  If  k=0  the  motion  is  independent  of  x,  i.e. 

{5}  =  a{v}  exp(iiot)  (6.30a) 

and  degenerates  to  a  one-dimensional  standing  wave.  This  wave  type  can 
occur  only  at  certain  frequencies  of  a  layered  medium  without  viscous 
damping.  The  natural  frequencies  can  be  found  from  the  secular  equation 

| [G]  -  u2[M] |  =0  (6.31) 

which  follows  from  Eqs.  (6.24)  and  (6.26).  When  ’"=0,  Eq.  (6.28)  shows 
that  the  horizontal  and  vertical  displacements  are  uncoupled.  The  motion 
corresponding  to  the  natural  frequencies  therefore  consists  either  of 
standing  S-waves  with  horizontal  displacements  or  of  standing  P-wave  with 
vertical  displacements.  The  standing  S-waves  are  identical  to  the 
degenerated  Love  waves  with  a  vanishing  wave  number  discussed  in  Chapter 
4,  except  that  the  displacements  occur  here  in  the  x-direction. 


Waves  of  the  types  b)  and  d)  do  not  exist  in  the  damped  case 
because  all  waves  attenuate  with  distance  from  their  source  due  to  damp¬ 
ing  and  thus  k2  cannot  vanish. 

A  study  of  the  structure  of  Eq.  (6.26)  with  reference  to  Eq.  (6.28) 
shows  that  if  k  is  an  eigenvalue  of  Eq.  (6.26)  and  (v)  the  corresponding 
eigenvector,  which  contains  the  elements  v  ,  j=l,  ...»  2n,  then  -k  is 
also  an  eigenvalue  and  its  corresponding  eigenvector  in  transposed  form 
is 


{v}  =  <-vx  v2  — v3  vA 


—v  V  > 

2n-l  2n 


(6.32) 


which  is  obtained  from  {v}  by  changing  the  sign  on  all  horizontal  dis¬ 
placements.  The  vector  {v}  is  the  adjoined  of  {v}. 

The  two  solutions  describe  the  motions 


{6}  =  a{v}  exp(iwt-ikx)  (6.33a) 

and 

{6}  =  a{v}  exp(iwt+ikx)  (6.33b) 
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which  are  identical  waves  propagating  or  decaying  in  opposite  directions. 
Substitution  of  the  second  solution,  -k  with  {v},  into  Eq.  (6.26)  yields 

([A]k2  -  i [B]k  +  [C]){vJ  =  {0}  (6.34) 

In  the  undamped  case,  when  [A],  [B]  and  [C]  are  real,  the  complex 
eigenvalues  occur  i.i  conjugate  pairs,  k  and  k.  For  this  case  the  complex 
conjugate  of  Eq.  (6.34)  is 

([A]k2  +  i [BJk  +  (CJ ) fv}  =  {0}  (6.35) 

which  shows  that  k  with  {v}  satisfies  Eq.  (6.26).  From  the  complex  con¬ 
jugate  of  Eq.  (6.26)  it  may  also  be  shown  that  -k  with  (v)  is  another 
solution. 

Summing  up,  the  eigenvalue  problem,  Eq .  (6.26),  has  the  property 
that  if 

a)  K  with  (v) 

is  a  solution  to  Eq.  (o.26)  then 

b)  -k  with  {v} 

is  another  solution  and,  in  the  undamped  case, 

c)  k  with  (  v) 

d)  -k  with  {v} 

are  also  solutions.  These  solutions  are  linearly  independent  with  the 
exception  that  the  solutions  c)  and  d)  depend  on  a)  and  b)  if  k  is 
purely  real  or  purely  imaginary. 


Orthogonality  of  Eigenvectors 

Two  solutions  to  Eq.  (6.26),  {v}^  with  k^O  and  {v}^  with  k^O,  are 
considered,  i.e. 

( [A]k2  +  i[B]k  +  [C])M  =  {0}  (6.36a) 

r  r  it 

and 

([A]k2  -  i[B]kg  +  [C]){v)s  =  {0}  (6.36b) 

Premultiplying  the  first  equation  by  (v)  /k^  and  the  second  by  {v}r/kg, 
then  subtracting  the  transposed  second  equation  from  the  first  gives 

{v}'J(([A!(kr-ks)  -  lC](kr-ks)/(krks)){v)r  =  0  (6.37) 
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Division  by  (k^-k^,  provided  k^kst  yields  the  orthogonality  relations 


{v}T[A]{v}  -  fv}T[C]fv}  /(k  k  )  = 

s  r  s  r  r  s 


2  if  k  =k 
r  s 

0  if  k  ?*k 
r  s 


(6.38) 


in  which  the  case  k^_-kg  is  used  to  normalize  the  eigenvectors.  Eq. 
(6.38)  may  be  expressed  by  che  matrix  equation 

[K] [V]T{A] [V] [K]  -  [V]T[C][V]  =  2[K]2  (6.39) 

in  which  [K]  is  a  2n  x  2n  diagonal  matrix  containing  one  each  of  the 
pairs  of  eigenvalues  ±kg,  s=l,  ...,  2n,  and  [V]  and  [V]  are  2n  x  2n 
modal  matrices  containing  the  corresponding  mode  shapes  (v)s  and  {v}^, 
s=l.  ...»  2n,  columnwise. 


6.3  Dynamic  Stiffness  Matrix  for  Layered  Region 

The  "energy  conditions"  discussed  in  Section  4.1  demand  that  only 
those  2n  generalized  Rayleigh  waves  which  decay  or  propagate  energy  in 
the  positive  x-direction  are  chosen  from  among  the  4n  solutions  to  Eq. 
(6.26).  Thus,  among  the  waves  possessing  complex  wave  numbers,  those 
waves  which  have  a  wave  number  with  a  negative  imaginary  part  are 
selected.  If  wave  numbers  are  real  the  selection  must  be  based  on  the 
direction  of  energy  propagation.  While  for  Love  waves  the  directions  in 
which  the  energy  and  the  "phase"  travel  are  always  identical,  it  some¬ 
times  occurs  that  the  "phase"  of  Rayleigh  waves  travels  backwards,  i.e. 
the  directions  of  the  phase  velocity  c=u)/k  and  of  the  energy  propagation 
are  opposed.  This  phenomenon  was  also  observed  by  Lysmer  [26]. 

It  will  be  shown  in  Section  6.4  that  the  energy  transmission  is 
positive  for  waves  with  a  positive  wave  number  k_  and  a  corresponding 

3 

mode  shape  (v)  if  the  complex  conjugate  of  the  mode  shape,  {v}  ,  is 
s  s 

equal  to  the  adjoined  mode  sh  .pe  {v}g  after  normalization  by  Eq.  (6.38), 

i.e.  {v}  =  {v}  .  If  {v}  =  -{v}  ,  the  direction  of  energy  transmission 

s  s  s  s 

is  opposite  to  that  of  the  phase  velocity.  Hence,  the  rules  for  the 
selection  of  the  2n  wave  numbers  ks>  s=l,  ...,  2n,  from  the  4n  wave 
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-  -  - 


s  s 

numbers  ±k  =  * (k  +  i  k~)  are 
S  i  z 


+  k  if  k  <  0  or  if  k0  =  0  and  {v}  =  {v} 

k  =  s  2  2  _  s  s 

-  k  if  k^  >  0  or  if  k*  -  0  and  {v}  =-{v) 

is/  /  s  s 


(6.40) 


The  displacements  in  the  layered  region,  when  expanded  into  waves 

corresponding  to  the  selected  wave  numbers  k  ,  are 

s 


{5}  =  ^  {v)s  exp(icot-iksx) 


(6.41) 


and  the  displacements  u.,  j=l,  ...,  2n,  of  the  nodal  points  at  the 
boundary  x=0  are,  omitting  the  time  factor  exp(iwt). 


(u}R  =  ^  {v}  a  =  [V]{al 

/  -<  s  s 

S=1 


(6.42) 


in  which  the  vectors  {u}  and  {a}  are  of  equal  dimension,  2n,  and  contain 
the  elements  u.  and  a^,  respectively.  The  matrix  [V]  is  the  2n  x  2n  modal 
matrix  containing  the  mode  shapes  (v)s  columnwise.  The  inversion  of 
Eq.  (6.42)  yields  the  mode  participation  factors  in  terms  of  the  nodal 
displacements,  i.e. 


la)  =  [v]-1{u}R 


(6.43) 


The  normal  stress  a  ,  and  the  shear  stress  t  at  the  boundary  x=0 
X  xz 

are  by  Eqs.(6.16)  and  (6.41)  for  z^<z<z^+^ 


°K  ■  1  [-iks<y2G.>1(2j+rz)/h.  •  v2j_1+<z"zj)/hj  •  Vy+l1 

s  =  l  , 

+  X./h.  •  (-v^.+v®.  .)  a 
3  J  2j  2j+2  1  s 


(fa. 44a) 


T*z*  SGj/hj  '  l(-V2j-l+v2j+l>  '  lksI(zj+r2>  *2j  +  <Z-Zj'V2i+21 1  % 


(6.44b) 
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in  which  vf  is  the  jth  element  of  {v}  .  The  time  factor  exp(iu)t)  is 
J  s 

implied. 

The  nodal  forces  which  are  in  equilibrium  with  these  stresses  and 
act  on  the  region  x>0  as  shown  in  Fig.  15b  are 


p  ’ 

2j-l 


s 

V2j-1 


2j  =  y  (i[A].k  +  fD],)  Vfj 

JII  w  J  s  J  s 

’  21+1  s=l  21 


P" 

2j+2 


V2j+1 

s 

V 

2j+2 . 


(6.45) 


in  which  [A],  and  [D]^  are  the  layer  matrices  defined  by  Eq.  (6.22).  The 
total  forces,  P  ^ ,  acting  on  the  layered  region  each  consist  of  the  force 
components  from  the  layer  above  and  below  the  respective  nodal  point. 


P  =  P1  +  P" 

j  j  j 


j — 1,  ....  2n 


with  Pj‘  =  P"  =  0.  Expressed  in  matrix  notation  they  are 
(P}R  =  (i f A] [V] [K]  +  [D] [V]){a> 


(6.46) 


(6.47) 


in  which  {p)K  is  the  vector  containing  the  elements  P^,  j=l,  ...,  2n, 
and  [K]  is  the  diagonal  matrix  containing  the  wave  numbers  kg,  s=l,  .. 


Substitution  of  Eq.  (6.43)  into  Eq.  (6.47)  gives 
{P}R=  [R] (u}R 


(6.48) 


in  which 


[RJ  =  i[A]  [V]  [K]  [ V] — A  +  [D] 


(6.49) 


The  2n  x  2n  matrix  [Rj  is  the  dynamic  stiffness  matrix  of  the  semi¬ 
infinite  layered  region  R.  It  relates  the  nodal  forces  to  the  simultaneous 
nodal  displacement  at  the  boundary  x=0. 

Matrix  [R]  is  symmetric,  as  it  should  be,  according  to  the  theorem 
of  reciprocity  [14 J.  Because  the  symmetry  is  not  apparent  from  Eq.  (6.49) 
it  is  proved  below. 

The  analysis  of  a  left  layered  region  I.  shown  in  Fig.  4  is  analogous 
to  that  of  a  right  layered  region  R.  The  only  difference  between  a  left 
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and  a  similar  right  layered  region  is  due  to  their  positions  with  respect 
to  the  global  coordinate  system.  Thus  the  dynamic  stiffness  matrix, 

[ L 1 ,  of  a  left  layered  region  may  be  computed  from  the  right  hand  side 
of  Eq.  (6.49)  by  changing  only  the  sign  on  all  the  coefficients  that 
relate  horizontal  forces  to  vertical  displacements  or  vertical  forces  to 
horizontal  displacements. 


Proof  that  [R]  is  Symmetric 

Eqs.  (6.49)  and  (6.25)  yield 

[R]T-[R]  =  i[V]~T[K][V]T[A]  +  [D)T  -  ifAHVHKHvr1  -  [D] 

=  i[V]“T[X] [V]-1  (6.50) 

in  which 

[xl  =  [K][V}T[A][V]  -  [V]T[A][V][K]  -  i[V]T[B)[V]  (6.51) 

In  order  to  show  that  each  element  x  of  [X]  vanishes,  two  solutions 

J.  b 

to  Eq.  (6.26)  are  considered,  i.e. 

([A]kg  +  i[B]kg  +  [C]){v)s  =  {0}  (6.52a) 

and 

( [A] k^  +  i[B]kr  +  [C]){v}r  =  {0}  (6.52b) 

T  T 

Pre-raultiplying  the  first  equation  by  tv}  and  the  second  by  {v  ,  then 
subtracting  the  transposed  second  equation  from  the  first  yields 

(v}'J  ( [A]  (kg-k^)  +  i[B](ks+kr)){v}s  =  0  (6.53) 

Since  (k  +kr)^0  for  any  two  k  of  matrix  [K],  Eq.  (6.53)  gives 

{v}^  ^[A](kr~ks)  -  i[B]jtv}g  =  0  (6.54) 

The  left  hand  side  of  Eq.  (6.54)  is  identical  to  element  xrg.  Hence, 
it  is  shown  that 

[R]  =  [R]T  (6.55) 
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"k  k 

Pre-mult iplying  the  first  equation  by  (v)  ,  the  second  by  {v}  ,  then 

r  s 

subtracting  the  conjugate  transposed  of  the  second  equation  from  the 
first  yields 


{v}*([A](kJ  -  kj)  +  i[B](kg  -  kr)){v}s  *  0  (6.62) 

Division  by  (k  -  k^), provided  k^k^,  proves  Eq.  (6.58). 

In  evaluating  Eq.  (6.60)  for  the  wave  types  discussed  in  Section 
6.2  it  is  of  assistance  to  consider  the  conditions  under  which  the  wave 
number  k  is  complex,  real  or  imaginary  in  the  undamped  case.  Pre- 

r  t  * 

multiplying  Eq.  (6.26)  with  lv)  and  introducing  the  abbreviations 

A  =  {v}  [A]{v)  ;  B  =  -i(v}  [B]{v}  ;  C  =  (v)  [C]{v)  (6.63) 

gives  the  quadratic  equation 


Ak  -  B  k  +  C  =  0 


which  has  the  solutions 


(b  ±  JT-  4AC  )  I  (2k) 


(6.64) 


(6.65) 


As  [A]  and  [C]  are  real  symmetric  and  [B]  is  real  skew  symmetric,  the 
numbers  A,  B  and  C  are  real,  A  being  positive,  because  matrix  [A]  is  positive 
definite.  With  the  above  abbreviations,  Eq.(6.60)  can  be  written 


Es  =  !Kl2(<ks  +  VA'B) 


(6.66) 


The  results  obtained  from  Eqs.  (6.65)  and  (6.66)  are: 


a)  k  is  complex  if  B^O  and  B~-4AC<0  and  from  k.  +  kg  -  B/A  it 


follows  that 


E  =  0 
s 


b)  k  is  imaginary  if  B=0  and  C>0  and  from  k  +  k  “  0  it  follows 


E  =  0 
s 


c)  k  is  real  if  B-4AC>0  and  from  (k  +  k  )  =  2k  and  Eq.  (6.66) 


it  follows  that 


to 

E  ~7 
s  4 


!  at  1 2(Ak  -  C/k  )  =  |  |<x  |2{v}*([A]k  -  [C]/k  ){v) 


(6.67) 
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In  Section  6. 2  it  was  shown  that  the  horitontal  and  vertical  dis- 

:i™  9o°  ~  *■  -  -  -  -  £ 

Eqs.  (6.67)  and 


case  fv}g-  ±{v>s  after  normalization  by  Eq.  (6.33) 

Oo\  .  .  -  / 


S  S 

(6.38)  therefore  yield 


E  = 
s 


r  it) 

+  2  ks  1%!  if  =  3-  {v} 

to  ,  ,  ,2 

~2ks  Kf  if  tv>s  =  -  {v} 


(6.68) 


which  is  of  the  same  form  as  the 
transmission  of  Love  waves,  Eq. 


corresponding  expression  for  the 
(4.58). 


energy 
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7.  AXISYMMETRIC  RAYLEIGH  WAVE  MOTION 


•  1  Motion  in  the  Layered  Region 


The  analysis  of  axisymmetric  Rayleigh  wave  motion  in  a  layered 
'  gion  of  the  type  shown  in  Fig.  5  is  similar  to  that  for  plane  Rayleigh 
wave  motion  presented  in  Chapter  6. 

Any  point  (r,z)  in  the  region  undergoes  displacements  in  the  radial 
direction,  r,  and  in  the  vertical  direction,  z.  For  harmonic  motion  at 
tne  frequency  <o  the  displacements  can  be  expressed  by 

Or(r,z,t)  =  qr(r,z,  exp(icot)  (7.1a) 

and 

Sz(r,z,t)  =  qz(r,z)  exp(icot)  (7.1b) 

The  governing  homogeneous  differential  equations  for  the  spatial 
part  of  the  displacements  are,  by  Eq.  (3.16),  for  the  jth  layer 


(>j  +  2Gj) 


and 


t2  Me)  +Cj\3z2  ~Mr) 

2 

+  p  .co  a  =0 
J  *r 


(7.2a) 


9r  i  ^q  3  q  \  /  S^q  ,  3q  3q  3q  \ 

j  j  \9z3r  r  oz  ^^2  j  ^  \  3r2  "*  r  "r  Zz^r  r  3z  J 


+  p  .co  q  =0 
•>  z 


(7.2b) 


Separation  of  the  variables  r  and  z  leads  to  a  particular  solu' ion 
to  Eqs.  (7.2)  of  the  form 


qr(r,z)  =  u(z)  H^2)(kr) 
q  ( r , z)  =  w(z)  i  l^2\kr) 

2  O 


(7.3a) 

(7.3b) 


(2)  (2) 

in  which  and  11^  are  the  Hankel  functions  which  were  introduced  in 
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Chapter  5.  The  functions  u(z;  and  w(z)  are  unknown  functions  which  must 
satisfy  the  coupled  ordinary  differential  equations  for  j=l,  ....  n 


2 

{ k2  ( X  .+2C .  1  -  u)2o.}u  -G.  -Mr  +  ik(X.+G.)  ^  =  0 
J  J  J  J  dz2  J  J  dz 


(7.4a) 


2 

{k2G.-a)2p  }  w  -  (A.+2G.)  ~  +  ik(A.-»G.)  ^  =  0 
J  J  J  J  dz2  J  J  ^ 


(7.4b) 


Eqs.  (7.3)  and  (7.4)  are  easily  verified  by  substituting  Eqs.  (7.3)  into 
Eqs.  (7.2)  and  observing  Eqs.  (5.8).  Eqs.  (7.4)  are  identical  to  the 
corresponding  equations  for  plane  Rayleigh  waves,  Eqs.  (6.4). 

The  boundary  conditions  at  the  horizontal  planes  z  ^ ,  j=l,  n+1. 


u(zj)  =  u(zp  ;  wfzp  =  wfzV) 
u(z  )  =  0  ;  w(?  )  =  0 


(7.5) 


a  (zV)  =0  ;  T  (z")  =  0 

z  1  zr  1 

o  (z!)  =  o  (z V)  ;  T  (z!)  =  T  (z ’.’) 
z  j  zv  y  zr  j  zr  j' 


When  the  stress  displacement  relations 


do  /  36  6  \ 

°z  =  (A  +  2G)  a r  +  x  + 

=  ((\  +  2G)  —  -  k  X  u  j  i  (kr)  exp(iiat)  (7.7a) 

/ 3  6  36  \ 

T  =  G  t  ~  +  rc  — ) 
zr  3z  ov 

=  G  (^  -  ikw)  H^(kr)  exp(itdt)  (7.7b) 


(7.6) 


are  substituted  into  L'q.  (7.6),  these  boundary  conditions  reduce  to  those 
of  the  plane  case,  Eq.  (6.9). 

As  the  differential  eqe  fions,  Eqs.  (6.4)  or  (7.4), and  the  boundary 
conditions,  Eq.  (6.9).  are  id*-  ical  for  plane  and  axisyraraetric  Rayleigh 
waves,  the  resulting  eigenvalue  problem  r  .d  its  solutions  arc  the 
same  in  both  the  plane  and  axisymmetric  case.  This  holds  not  only  for 
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the  continuum  theory  but  also  for  the  discrete  theory,  because  the 
layered  region  is  discretized  in  the  axisymmetric  case  in  the  same  manner 
as  in  the  plane  case. 

It  is  again  assumed  that  u(z)  and  w(z)  vary  linearly  within  each 
layer  according  to  Eq.  (6.14).  The  displacements  in  the  layerec  region 
can  then  be  defined  by  the  displacements  of  the  layer  interface.;  which 
are  collected  in  the  vector  {5},  i.e. 


{6}T  =  <5  <5  6  „  5  .  . . .  6  6  > 

rl  zl  r2  z2  rn  zn 


(7.8) 


in  which  6  .  and  6  .  are  the  radial  and  the  vertical  displacements  of 
rj  zj 

the  interfaces. 

The  displacements  may  be  expressed  as  a  linear  combination  of  the 
2n  modes  that  propagate  or  decay  in  the  positive  r-direction,  i.e. 


{5}  =  £  {6 } 


(7.9) 


in  which  is  the  participation  factor  of  the  sth  mode  {6}s. 

The  sth  mode  {6}g  is  obtained  from  the  sth  solution  to  Eq.  (6.26) 
i.e.  from  (v)s  and  ks>  when  Eqs.  (7.3)  are  observed.  The  displacements 
of  the  sth  mode  are  accordingly 


0s  =  H^(k  r)  •  exp(iwt) 
J  2  1  s 


6s  =  v®iH^(k  r)  •  exp(iiot) 


j=l,3,  ...»  2n-l  (7.10a) 


j=2,4,  ...,  2n 


(7.10b) 


The  2n  modes  used  in  the  expansion,  Eq.  (7.9), are  selected  accord¬ 
ing  to  the  rules  of  Eq.  (6.40)  which  were  derived  for  the  plane  case. 
From  the  properties  of  tue  Hankel  functions  and  from  the  similarity  of 
axisymmetric  waves  to  plane  waves  at  large  values  of  |kgr|  it  follows 
that  these  rules  also  hold  for  the  axisvmmetric  case. 
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and 


i  =  G.  l-j-  u  (z)  -  ik  w  (z)\  a  a 
rz  Z-  j\dz  s  s  s  /  ss 


s=l 


in  which 


u  ( z)  =  (z.,_  -  z)/h.  •  v„.  ,  +  (z-z.)/h.  •  V- . 


'j+l 


j  2j-l 


and 


ws(z)  =  (zj+1  -  z)/h.  •  v2j  +  (z-zp/h^  •  v0. 


(7.16b) 

s 

2  j+l 

(7.17a) 

s 

2  j+2 

(7.17b) 

The  nodal  forces  which  are  in  equilibrium  with  the  stresses  on  a 
segment  of  unit  raaian  are 


P' 

2j-l 

s 

v  • 

2j-l 

b 

s 

^2j-l 

n  • 

2n 

S 

2n 

,S 

23 

P2j+1 

•  =  %  ik  r  [A] .  ‘ 
ZL  so  j 

S=1 

V2j 

vS  * 

2j+l 

a 

s 

b 

s 

'  a  +  y  r  ([D]  +[EJ  )  < 

25  A-—/  U  J  J 

S  =  1 

♦;j+i 

.  P2j+2 

vS  • 

2  j+2 

a 

s 

.  *^2  j+2  . 

>  a 


(7.18) 

in  which  the  layer  matrices  [A],  and  [D] .  are  those  defined  by  Eqs.  (6.22) 
and  layer  matrix  [E]^  is 


i  G- 

[E ] .  = 

J  3  ro 


2  0  10 
0  0  0  0 
10  2  0 
0  0  0  0 


(7.19) 


The  total  forces  acting  on  the  layered  region  r>ro  at  the  boundary 
nodal  joints  each  consist  of  the  force  contributions  from  the  layer 
above  and  below  the  respective  nodal  joint,  i.e. 

(7.20) 


p  =  p»  +  p" 

.1  j  j 


^“1,  ..*,  2n 


with  P’t  =  P"  =  0.  The  forces  expressed  in  matrix  notation  are 

{P }R  =  ro(i[A][Y][K]  +  [DU*]  +  [EjK>|‘:  i)  (7.21) 
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in  which  the  vector  {p}  contains  the  nodal  forces  P  ,  j=l,  . ...  2n. 
Matrices  [A],  [K]  and  [D]  were  defined  in  Chapter  6.  Matrix  [E]  is 
assembled  by  addition  of  the  layer  matrices  fEj .  as  shown  in  Fig.  14 
by  substituting  [Ej^  for  [Xj^-  Matrix  [4']  consist  of  the  column  vectors 
{^}s>  s=l ,  ....  2n,  which  have  the  elements 


,s  s  , 
>b .  =  v .  b 


j=l,  3,  . ..,  2n-l 


(7.22a) 


<pS  =  vS  a  j=2,  4,  ....  2n  (7.22b) 

J  J  s 

Substitution  of  Eq.  (7.14)  into  Eq.  (7.21)  gives 

{p)R=[R]{u}R  (7.23) 


in  which  [R]  is  the  dynamic  stiffness  matrix  for  the  semi-infinite 
layered  region  R,  i.e. 

[R]  =  r^ilAjmiKH*]"1  +  [DJ  +  [EJ)  (7.24) 

which  relates  the  nodal  forces  per  radian  to  the  simultaneous  nodal  dis¬ 
placements  at  the  boundary  r3^.  The  matrix  [R]  has  properties  similar 
to  those  of  the  dynamic  stiffness  matrix  for  the  plane  case;  in  parti¬ 
cular,  matrix  [R]  is  symmetric  in  accordance  with  the  theorem  of 
reciprocity  [14].  This  may  be  shown  by  a  proof  similar  to  that  presented 
for  the  plane  case  in  Section  6.3. 
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8.  SUMMARY  OF  PROCEDURE 


The  theory  for  the  proposed  method  of  analysis  is  now  completed 
and  the  steps  involved  in  analyzing  plane  or  axisymmetric  systems  of 
the  type  shown  in  Figs.  4  and  5  are  summarized  below. 

i .  idealization 

1.  Subdivide  the  layered  regions  R  ar.d  L  into  sublayers  so  that 
the  thickness  of  each  sublayer  is  less  than  about  1/10  of 
the  length  of  S-waves  ^r; veiling  in  it  at  the  frequencies 
considered. 

2.  Cover  the  irregular  region  I  by  a  finite  element  mesh  and 
specify  the  coordinates  and  reference  numbers  of  the  nodal 
joints.  The  nodal  joints  at  the  boundaries  between  the 
irregular  region  and  the  layered  regions  must  coincide  with 
the  interfaces  of  the  sublayers,  see  Fig.  6.  The  maximum 
dimension  of  each  element  should  again  be  less  than  1/10  of 
the  length  of  S-waves  and  perhaps  smaller  near  stress 
concentrations . 

3.  Define  the  material  properties  of  the  finite  elements  and 
sublayers.  Viscous  damping  is  included  by  the  use  of  complex 
moduli. 

4.  Define  the  prescribed  nodal  displacements  and  the  external 
forces  acting  on  nodal  joints  of  Region  I  and  specify  the 
frequency  of  the  harmonic  loading. 


j i.  Dynamic  Stiffness  Matrix  for  Layered  Regions 

1.  Formulate  the  eigenvalue  problem  for  region  R,  i.e.  Eis.  (4.31) 
and  (6.26)  in  the  Love  and  Rayleigh  wave  cases,  respectively. 

2.  Solve  the  eigenvalue  problem  by  the  methods  outlined  in 
Appendix  1  and  2. 
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3.  Form  the  dynamic  stiffness  matrix  [R]  for  region  R  according 
to  Eqs.  (4.47).  (5.28),  (6.49)  and  (7.24)  in  the  plane  and 
axisymmetric  Love  and  Rayleigh  wave  cases,  respectively. 

4.  In  the  plane  cases,  perform  steps  1.,  2.  and  3.  in  a  similar 
manner  to  obtain  Lhe  dynamic  stiffness  matrix  [L]  for  region  L. 

Ill .  Analysis  of  Irregular  Region 

1.  Form  element  mass  and  stiffness  matrices  [M']  and  [K'J.  Use 
complex  moduli  to  include  viscous  damping  in  [K*]. 

2.  Assemble  the  global  matrix  [A]  =  [K]  -  to"  [M]  by  the  direct 
stiffness  method  from  the  element  matrices  of  region  I  and  the 
frequency  dependent  dynamic  stiffness  matrices  [R]  and  [L]  for 
the  layered  regions. 

3.  Compute  the  nodal  displacements  by  solving  Eq.  (3.27).  The 
amplitudes  and  phase  angles  of  the  nodal  displacements  are 
determined  by  Eqs.  (3.5)  and  (3.6).  Reaction  forces  may  be 
computed  by  Eq.  (3.31). 


IV.  Displacements  in  Layered  Regions 

1.  Compute  mode  participation  factors  for  region  R  by  Eqs.  (4.37), 
(5.19),  (6-43)  and  (7.14)  in  the  plane  axisymmetric  Love 
and  Rayleigh  wave  cases,  respectively. 

2  Compute  displacements  of  the  sublayer  interfaces  in  region  K 
by  Eqs.  (4.41),  (5.17),  (6.41)  or  (7.9). 

3.  In  the  plane  cases,  perform  steps  ].  and  2.  for  region  L. 


Computer  programs  which  accomplish  steps  II,  III  and  IV  of  the 
above  procedure  for  systems  c v  the  type  shown  in  Figs.  4  and  5  have  been 
developed.  The  major  part  of  the  computation  required  is  performed  in 
complex  variables  utilizing  the  complex  arithmetic  capabilities  of 
FORTRAN  IV.  Res>’lts  from  the  analysis  of  several  examples  are  presented 
in  the  following  chapters. 
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9 .  LINE  LOADS  ON  HOMOGENEOUS  LAYER 


The  accuracy  of  Che  presented  method  is  studied  by  application  to 
the  following  problem.  A  harmonic  line  load  of  unit  amplitude  acts  on 
the  surface  of  a  homogeneous  layer  of  unit  thickness  as  shown  in  Fig.  16. 
The  layer  material  has  a  unit  mass  density  and  a  Poisson's  ratio  v=0.3. 

I f  no  material  damping  is  assumed,  the  shear  modulus  G=1  unit,  and  if  a 
fraction  of  critica.  damping  $=0.05  is  chosen,  the  stress-strain  behavior 
is  represented  by  the  complex  shear  modulus  G=(l  +  i2£)  units. 


9. 1  Love  Wave  Case 


Fig.  16  shows  <i  line  load,  P=1  *exp(i!j*:)  per  unit  length,  which 
acts  in  the  direction  perpendicular  to  the  specified  x-z-plane  and  thus 
generates  plant  generalized  Love  waves.  The  correct  solutioi  to  this 
problem,  labelled  E,  is  derived  in  Appendix  3  and  serves  to  check  the 
numerical  solutions  which  are  computed  by  the  presented  method  for 
three  different  discrete  models.  Each  of  these  models  makes  use  of  the 
symmetry  of  the  problem  about  the  z-axis. 

Model  A,  shown  in  Fig.  17,  consists  of  a  region  I  with  18  x  18 
rectangular  finite  elements  and  a  region  R  with  L8  sublayers.  The 
finite  element  mesh  is  condensed  near  the  line  load  and  the  maximum  mesh 
size  is  1/15  of  the  wave  length  of  shear  waves  at  the  highest  frequency 
considered  (u-2 “). 

Model  B  contains  no  finite  element  region;  the  layered  region  R 
extends  from  x=Q  to  infinity  and  consists  of  18  sublayers  as  does 
region  R  in  Model  A. 

Model  C  is  obtained  from  Model  B  by  dividing  each  sublayer  into 
two  sublayers  of  equal  thickness  resulting  in  a  total  of  36  sublayers. 

The  solutions  ootained  from  these  models  are  referred  to  as 
solutions  A,  B,  and  C  respectively. 

Displacements  for  the  static  case,  w=0,  ana  displacement  amplitudes 
for  harmonic  motion  with  and  without  material  damping  at  the  frequency 
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to=2ir  are  shown  in  Table  1  for  several  points.  The  phase  shifts  of  the 
displacements  with  respect  to  the  phase  of  the  line  load  are  presented 
in  Table  2. 

Table  1  and  Table  2  show  that  the  solutions  A  and  B  are  about 
equally  close  to  the  correct  solutions,  E,  and  satisfactorily  accurate 
in  the  neighborhood  of  the  line  load  and  at  some  distance  from  it.  The 
solution  C  is  significantly  better  than  the  solution  B  because  the  sub¬ 
layers  in  Model  C  are  only  half  as  thick  as  in  Model  B. 

The  errors  in  the  displacements,  displacement  amplitudes  and  phase 

shifts  of  the  solution  C  are  almost  precisely  four  times  smaller  than 

the  errors  of  the  solution  B.  This  is  true  for  all  the  values  presented 

in  Table  2  and  Table  2  and  is  therefore  a  strong  indication  that  the 

2 

errors,  if  they  are  small,  are  essentially  proportional  to  1/h  where  h 
is  the  thickness  of  a  sublayer. 

This  observation  can  be  utilized  to  obtain  a  much  improved  solution, 
called  D,  from  the  solu  ions  B  and  C  by  adding  one  third  of  the 
difference  between  the  solutions  C  and  B  to  the  solution  C.  Values  of 
the  solution  D  are  also  listed  in  Table  1  and  Table  2.  They  have  one  or 
two  more  correct  digits  than  the  values  of  solution  C  and  agree  very 
well  with  the  values  of  the  correct  solution. 

The  approximate  solutions  for  the  layered  regions  consist  of 
finite  series  of  approximate  wave  modes,  Eq.  (4. 41),  while  the  correct 
solutions  are  infinite  series  of  exa^t  "ave  modes,  Eq.  (A3. 13).  Some 
wave  numbers  of  corresponding  approximate  and  exact  modes  are  listed  in 
Tabic  3,  in  which  the  letters  B,  C,  D,  and  E  indicate,  as  before,  the 
respective  models  an’  solutions. 

The  values  C,  obtained  with  36  layers,  are  again  almost  exactly 
four  times  closer  to  the  correct  values  E  than  the  values  B  obtained 
with  18  layers.  The  very  precise  values  D  are  computed  by  adding  one 
third  of  the  difference  between  the  values  C  and  B  to  the  values  C. 

One  may  recall  that  the  displacements  of  the  sth  wave  mode  are 
6s(x,z,t)  =  vg(z)  exp(iwt-iksx) .  In  the  undamped  case,  3=0,  the  wave 
numbers  kg  are  either  real  or  imaginary  and  correspond  to  waves  which 
propagate  with  constant  amplitude  or  decay  exponentially  in  the  x- 
direction.  The  rapidly  decaying  waves,  which  have  a  wave  number  with 
a  large  negative  imaginary  part,  contribute  significantly  to  the  total 
displacements  only  at  small  values  of  x,  i.e.  close  to  the  line  load. 
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At  some  distance  from  the  line  load  only  the  propagating  modes  are 
important. 

When  damping  is  irtroduced,  as  in  the  case  B=0.05,  all  waves 
propagate  and  attenuate  in  the  x-direction  because  all  wave  numbers  are 
complex.  However,  the  exponential  decay  is  relatively  small  for  the 
waves  which  would  propagate  with  constant  amplitude  in  the  undamped 
case  and  the  wave  lengths  of  the  modes  which  wouLd  decay  exponentially 
in  the  undamped  case  are  relatively  long,  as  the  real  parts  of  their 
wave  numbers  are  small. 


Rayleigh  Wave  Case 


The  case  in  which  the  harmonic  line  load  P=l*  exp  (lent)  per  unit 
length  is  applied  at  x=z=0  and  acts  in  the  z-direction,  rather  than  in 
the  y-direction  as  shown  in  Fig.  16,  is  now  considered.  Generalized 
plane  Rayleigh  waves  are  generated  by  the  force. 

The  correct  solution  to  this  problem  is  not  known  and  cannot  be 
found  easily.  Approximate  solutions  are  obtained  for  the  discrete 
models  A,  B,  and  C  described  in  the  Love  wave  case.  Their  values  for 
displacements,  displacement  amplitudes  and  phase  shifts  are  shown  in 
Table  5  and  Table  6  for  several  points. 

The  agreement  between  the  solutions  A,  B,  and  C  is  about  as  good 
as  in  the  Love  wave  case.  It  shows  that  the  numerical  results  change  only 
very  slightly  as  the  left  boundary  of  the  layered  region  R  is  moved  from 
x=l  (Model  A)  to  x=0  (Mode.1  B)  and  as  the  discretization  of  the  region 
R  is  refined  by  reducing  the  thickness  of  the  sublayers  (Model  C) .  The 
solutions  dxffer  significantly  oii’y  at  points  very  close  to  the  line 
load,  see  point  (x=0.01,  z=0) .  Tais  must  be  expected,  because  the 
strain  directly  under  the  line  load  (x=r=0)  is  infinitivelv  large  in  a 
linearly  elastic  continuum. 

Some  approximate  wave  numbers  which  were  computed  from  Model  B 
and  C  are  compared  with  correct  wave  numbers  in  Table  4.  The  correct 
numbers  are  the  roots  of  a  difficult  secular  equation  with  transcendental 
terras  [22].  They  were  computed  iteratively  by  the  secant  method  using 
as  initial  guess  the  wave  numbers  obtained  from  Model  C. 

The  values  C,  obtained  with  36  layers,  are  again  almost  exactly 
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four  times  closer  to  the  correct  values  E  than  are  the  values  B,  obtained 

with  18  layers.  As  in  the  Love  wave  case,  this  indicates  that  the  errors, 

2 

ii  they  are  small,  are  essentially  proportional  to  1/h  where  h  is  the 
thickness  of  a  sublayer. 

The  values  D  are  computed  by  adding  one  third  of  the  difference 
between  the  values  C  and  B  to  the  values  C.  This  extrapolation  improves 
the  accuracy  by  one  or  two  digits. 
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10.  SCREENING  EFFECT  OF  TRENCHES 

Many  attempts  have  been  made  to  reduce  ground-transmitted  vibrations 
by  the  installation  of  wave  barri**  -s  in  the  form  of  trenches  or  walls. 
Barkan  [9]  reviewed  several  installations  and  came  to  the  conclusion 
that  such  barriers  are  often  practically  useless,  ne  attributed  the 
iimited  success  to  the  absence  of  a  rational  design  procedure  for  this 
tyue  of  installation.  Recent  experimental  studies  by  Woods  [52]  and  by 
Dolling  [15]  have  greatly  improved  the  understanding  of  the  screening 
problem.  Bv.  a  rational  design  procedure  which  includes  a  sound  analysis 
of  the  screening  is  still  not  available. 

Tne  method  developed  in  this  dissertation  is  capable  of  analyzing 
plane  strain  and  axisymmetric  screening  problems.  A  relatively  simple 
example  is  presented  below  in  which  the  ground  transmitted  disturbance 
consists  of  plane  generalized  Love  waves. 

Figs.  18  and  19  show  computed  surface  displacements  near  a  st^ip 

footing  embedded  in  a  homogeneous  linearly  elastic  soil  layer  over  a 

2 

rigid  base.  The  footing  is  excited  by  a  harmonic  force  P  =  w  •  coscot 
per  unit  length.  This  type  of  loading  is  typical  for  foundations  support¬ 
ing  rotatory  machines. 

The  finite  element  mesh  used  in  the  a  alvsis  contained  450 
rectangular  elements  and  extended  from  the  left  edge  or  the  footing  to 
the  right  edge  of  the  trench.  The  semi-infinite  layered  regions  c>  the 
left  and  to  the  right  of  the  finite  element  mesh  were  represented  by 
18  sublayers  each. 

Fig.  19  shows  the  effect  of  trench  depth  at  the  frequency  w  =  1. 

At  this  frequency  a  trench  is  quite  effective  in  reducing  the  surface 
vibrations.  Its  effectiveness  increases  with  increasing  trench  depth. 

A,  perhaps,  unexpected  result  is  that  the  presence  of  the  trench  causes 
a  reduction  in  the  displacement  amplitude  of  the  footing.  This 
reduction  is  due  to  reflections  from  the  trench  wall.  However,  the 
reflections  might  just  as  well  cause  an  amplification  of  the  footing 
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displacement  depending  on  the  frequency  of  the  exciting  force  and  on  the 
distance  between  the  footing  and  the  trench.  The  curve  labeled  a)  =  1.5 
in  Fig.  18  illustrates  this  point. 

Fig.  18  shows  the  effect  of  changing  the  frequency  for  the  case  in 
which  the  trench  penetrates  half  way  through  the  soil  layer.  At  the 
lower  frequencies  to  =  0.5  and  to  =  1.0  the  trench  is  quite  effective  but 
as  the  frequency  is  increased  to  to  =  1.5  the  trench  becomes  useless. 

This  phenomenon  is  associated  with  the  number  of  propagating  Love  modes 
in  the  system.  At  frequencies  below  to  =  ir/8  =.39  no  propagating  modes 
exist  and  the  displacement  amplitudes  in  the  far  field,  i.e.  beyond  the 
trench,  decay  exponentially  according  to  Eq.  (7;.38b).  In  the  frequency 
range  0.39  to  1.18  one  propagating  mode  exists.  Its  mode  shape,  shown 
in  Fig.  18,  indicates  that  most  of  the  energy  is  transmitted  in  the  upper 
half  of  the  layer.  The  trench  is  therefore  quite  effective  in  reflecting 
the  Love  wave.  This  is  demonstrated  by  the  curves  labeled  w  =  0.5  and 
1.0  in  Fig.  18.  As  the  frequency  is  increased  into  the  range  1.18  to 
1.96  a  second  propagating  mode  appears.  Its  mode  shape,  shown  in  ^ig.  18, 
indicates  that  a  significant  part  of  the  energy  is  transmitted  below 
the  bottom  of  the  trench.  Consequently,  the  trench  is  not  an  effective 
barrier  for  the  second  mode.  This  is  illustrated  by  the  curves  labeled 
(0  =  1.5  in  Fig.  18.  At  frequencies  beyond  u>  =  1.96  additional  modes 
will  appear  and  the  displacement  in  the  far  field  becomes  more  complicated. 
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11.  TORSIONAL  MOTION  OF  RIGID  CIRCULAR  FOOTINGS 


11.1  Introduction 

A  rigid  circular  footing  on  an  axisymmetric  medium  has  four 
distinct  degrees  of  freedom;  they  are  the  vertical  and  horizontal  trans¬ 
lations  and  the  rotations  about  a  horizontal  and  a  vertical  axis.  The 
rotational  motion  about  the  vertical  axis  of  symmetry,  called  torsional 
motion,  is  of  considerable  importance  to  precision  tracking  radar  towers 
and  communication  antennas  [47,  48]. 

Radar  towers  are  often  supported  on  spread  footings  of  approximately 
circular  plan.  The  footings  are  st:‘ff  compared  to  the  soil  underneath 
and  may  therefore  be  treated  as  rigid  circular  disks. 

A  solution  for  the  harmonic  torsional  motion  of  a  rigid  circular 
disk  fixed  to  the  surface  of  a  homogeneous  linearly  elastic  half  space. 
Fig.  1,  was  first  obtained  by  Reissner  and  Sagoci  [41] .  Stallybras 
[«4]  derived  a  solution  to  the  same  problem  which  is  easier  to  evaluate 
for  high  frequencies.  Weissmann  [47]  considered  the  effect  of  material 
damping  and  slippage  between  the  footing  and  the  surface  of  the  half 
space. 

The  assumption  that  the  footing  is  supported  by  a  homogeneous 
half  space  is  mathematically  convenient  and  for  practical  purposes 
suitable  if  the  subsoil  is  homogeneous  to  a  depth  of  a  few  footing 
diameters.  But  this  condition  is  not  often  encountered  in  practice. 

Arnold  et  al.  [6]  obtained  a  solution  for  the  torsional  motion  of 
a  rigid  circular  footing  supported  by  an  elastic  layer  which  extends 
laterally  to  infinity  and  is  fixed  at  its  bottom  to  a  rigid  base,  see 
Fig.  2.  They  assumed  that  the  dynamic  stress  distribution  under  the 
footing  on  a  layer  is  similar  to  the  static  stress  distribution  under  a 
rigid  footing  on  a  half  space.  This  assumption  loses  its  validity  as 
the  frequency  of  the  motion  increases  and  the  thickness  of  the  layer 
decreases.  Awojobi  [  8]  treated  the  same  problem  by  reducing  it  to  the 
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solution  of  an  integral  equation  without  assuming  a  priori  a  stress 
distribution  under  the  rigid  footing.  However,  he  substituted  the  kernel 
of  the  integral  by  an  approximation  for  which  he  could  find  a  closed  form 
solution. 

The  half  space  solution  based  on  the  model  in  Fig.  1  and  the  layer 
solution  based  on  the  model  in  Fig.  2  are  currently  used  in  design 
practice.  But  as  footings  are  always  embedded  in  the  ground  and  are 
usually  supported  by  layered  soil,  the  applicability  ui  the  half  space 
and  layer  solu^'ans  is  often  questionable.  No  methods  of  analysis  have 
been  published  so  far  which  can  quantitatively  account  for  the  effect  of 
embedment  and  layering  of  the  soil. 

the  method  presented  herein  has  the  flexibility  tc  analyze  the 
response  of  footings  embedded  in  layered  soil.  In  the  following  examples 
the  effect  of  the  "hickness  of  a  homogeneous  layer  over  a  rigid  base,  the 
effect  of  damping,  the  effect  of  embclr..eit  and  the  effect  of  an  increase 
in  the  shear  modulus  or  the  soil  with  •  uptb  are  studied. 

First,  a  few  terms  used  in  the  presentation  and  the  discussion  of 
the  examples  have  to  be  explained. 

The  harmonic  torsional  motion  of  a  massless  rigid  footing  is  des¬ 
cribed  by 

T  * 

9  =  •—  F  exp(iu>t)  (11.1s) 

e 

in  which  0  is  the  angle  of  rotation  about  the  axis,  Tq  is  the  amplitude 
of  the  external  torque  on  the  footing,  and  kg  is  the  static  spring 
constant.  The  dimensionless  complex  quantity  F  =  F^  +  i  Fj,  designated 
as  the  displacement  function,  depends  on  the  circular  frequency  w  or  on 
the  dimensionless  frequency  ratio 

a  =  cor  /V  (11.2) 

o  os 

in  which  rQ  is  the  radius  of  the  footing  and  V  is  the  shear  wave  velocity 
in  the  soil. 

If  the  displacement  function  F  for  a  footing  without  inertia  is 
known  one  can  calculate  the  response  of  a  footing  with  ’"he  mass  moment 
of  inertia, Iq, about  the  axis  of  rotation,  as  shown  by  Lysmer  and  Ricr.art 
[24],  by 
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0  =  r—  M  cos(wt  +  <)>) 


in  which  K  <s  the  dynamic  magnification  factor 


M  = 


2  2 

F  4-  F 
*1  2 


(1  -  </lQ  Fj/kg)*  +  (u  I0  F2/kQ)^ 
and  4>  is  the  phase  shift 


4-  =  tan 


-1 


h  -  “%  <Fi  +  F2)/I:e 


(11.3) 


(11.4) 


(11.5) 


It  is  often  convenient  to  employ  a  dimensionless  inertia  ratio, 
which  can  be  defined  by  142] 


Be 


(11.6) 


where  p  is  the  mass  density  of  the  soil. 

If  the  mass  density  and  the  shear  wave  velocity  of  the  soil  vary 
with  depth,  p  and  Vg  in  Eqs.  (11.2)  and  (11.6)  are  the  values  of  the 
density  and  the  shear  wave  velocity  at  a  designated  reference  depth.  In 
all  the  following  examples,  in  which  soil  properties  vary  with  depth, 
this  reference  depth  will  be  taken  to  be  equal  to  one  footing  radius. 


11.2  Comparison  of  Numerical  Results  with  Analytical  Solutions  and 
Experimental  Results 

Arnold  et  al.  [6]  performed  model  tests  to  check  their  analytical 
solutions  for  the  torsional  motion  of  a  rigid  circular  disk  on  an  elastic 
layer  against  experimental  results.  They  glued  a  disk  to  the  surface  of 
a  sheet  of  foam  rubber  and  forced  the  disk  to  vibrate  in  the  torsional 
mode  . 


5 


Fig.  20  shows  their  experimental  and  their  analytical  results  for 

the  following  parameters.  The  ratio  of  the  layer  thickness  to  the  disk 

radius  is  H/r  =  0.97  ar^  the  inertia  ratio  of  the  disk  is  Bft  =  0.79. 
o  u 
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The  analytical  and  experimental  resul ts  agree  well  with  respect  to  the 
resonance  frequency  and  the  general  shape  of  the  response  curve;  but 
they  differ  in  the  maximum  amplitude. 

Awojobi  f S ]  also  compared  his  analytical  solution  with  Arnold's 
experimental  results.  The  comparison,  see  Fig.  20,  shows  a  good  agree¬ 
ment  in  the  ascending  part  of  the  response  curve,  but  a  difference  of 
about  10  percent  in  the  resonance  frequency. 

The  same  example  was  studied  numerically  by  the  finite  element 
method  with  a  fine  mash,  which  is  shown  in  Fig.  21.  The  maximum  dimen¬ 
sion  of  the  elements  and  the  maximum  thickness  of  the  sublayers  are 
less  than  1/25  of  the  length  of  shear  waves  at  the  frequency  ratio 
aQ  =  2.5.  Numerical  results  are  presented  in  Fig.  20  for  an  elastic 
material  without  damping,  8=0%,  and  for  a  viscoelastic  material  with 
a  fraction  of  critical  damping  8  =  7.5%.  The  numerical  solution  for  the 
undamped  material  and  Arnold's  experimental  results  agree  very  well 
with  respect  to  the  resonance  frequency  and  the  general  shape  of  the 
response  curve.  The  deviation  in  the  peak  amplitude  is  apparently 
caused  by  the  internal  damping  of  the  foam  rubber.  The  numerical  solution 
for  8  =  7.5%,  which  takes  the  damping  of  the  foam  rubber  into  account, 
fits  the  experimental  results  exceptionally  well  throughout  the  given 
frequency  range. 

The  small  difference  between  Arnold's  analytical  solution  and  the 
finite  element  solution  for  8  =  0%  can  be  attributed  mainly  to  * 
deviation  of  the  stress  distribution  assumed  by  Arnold  et  al .  from  the 
actual  dynamic  stress  distribution  under  the  rigid  footing.  The  error 
resulting  from  the  finite  element  discretization  is  expected  oe  small 
because  the  chosen  mesh  is  fine. 

The  computational  time  needed  for  the  finite  elerar-.  malysis  at 
one  frequency  was  approximately  4  seconds  on  the  CD C  A.  computer. 


11.3  Effect  of  Layer  Thickness 

Fig.  22  shows  the  displacement  function  F  for  the  torsional  motion 
of  a  massless  footing  on  a  homogeneous  elastic  layer  without  material 


damping.  The  figure  also  shows  the  geometry  of  the  vibrating  system. 
The  radius  of  the  footing  is  rQ  and  the  thickness  of  the  layer  is  H. 
Displacement  functions  are  presented  for  the  ratios  H/rQ  =  1,  2,  4,  8, 
and  00 .  The  function  for  H/ =  00  is  the  Reissner-Sagoci  solution  for 
the  elastic  half  space  {41]. 

Fig.  22  shows  that  the  F2~curves  are  zero  for  frequency  ratios 
below  certain  values  which  can  be  computed  from  Eq.  (A3. 6).  No  Love 
waves  can  exist  which  propagate  through  the  layer  and  no  energy  is 
dissipated  away  from  the  footing,  if 


it  ro 
ao  2  H 


(11.7) 


l  =  27T  V  /oj  =  2ir  r  /a  >  4H 
s  s  o  o 


(11.8) 


in  which  Zs  is  the  length  of  a  shear  wave  and  Vg  is  the  shear  wave 
velocity.  This  means  that  there  is  no  geometric  or  radiation  damping, 
that  the  rotation  of  the  footing  is  in  phase  with  the  applied  torque, 
and  thaf  the  F~-curve  is  2ero.  In  this  range  where  Z  >  4H,  the  response 
of  a  footing  with  inertia  can  become  infinitely  large  as  the  denominator 
in  Eq.  (11.4)  may  vanish,  depending  upon  the  moment  of  inertia  of  the 
footing  and  on  the  frequency. 

The  waviness  of  the  displacement  function  Fig.  22  is  caused  by  the 
change  of  wave  modes  from  the  decaying  type  to  the  propagating  type  as 
the  frequency  ratio  increases,  see  Section  4.3.  The  change  of  the  sth 
mode  occurs  at 


2  n  H 


(11.9) 


according  to  Eq.  (A3. 6).  The  wave  number  kg  of  the  sth  mode  vanishes  at 
this  aQ  value,  while  it  is  purely  imaginary  and  purely  real  for  smaller 
and  for  larger  frequency  ratios,  respectively. 

As  pointed  out  in  Section  5.2,  the  dynamic  stiffness  of  a  layer 
against  torsional  motion  does  not  become  singular  if  one  of  the  wave 
numbers  vanishes. 

The  displacement  functions  for  small  ratios  H/rQ  differ  significantly 
from  the  half  space  solution;  the  smaller  H/rQ  the  greater  they  differ. 
However,  the  displacement  function  for  H/r^  =  8  shows  good  agreement  with 
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the  half  space  solution  in  the  range  aQ  <  2,  which  is  the  range  of  interest 
for  most  footings. 

The  values  of  the  static  spring  constants  are  also  presented  in 
Fig.  22.  For  H/rQ  >  2  the  static  spring  constant  is  close  to  that  of  a 
footing  on  a  half  space. 

The  discrete  models  used  in  the  computation  of  the  displacement 
functions  were  similar  to  the  one  shown  in  Fig.  21.  The  maximum  mesh 
size  was,  in  all  cases,  smaller  than  1/12  of  the  length  of  a  shear  wave. 
Values  of  the  displacement  functions  were  computed  at  frequency  ratios  in 
intervals  of  0.2. 

11.4  Effect  of  Embedment 

It  was  shown  in  the  previous  example  that  a  footing  on  a  homo¬ 
geneous  elastic  layer  with  a  thickness  of  several  footing  radii  responds 
to  torsional  excitation  very  much  the  same  as  does  a  footing  on  a  homo¬ 
geneous  elastic  half  space. 

The  effect  of  embedment  will  be  demonstrated  for  a  footing  on  a 
thick  layer,  H/r^  =  8.  The  response  curves  for  three  different  footings, 
each  with  inertia  ratios  Bg  =  2  and  Bg  =  4,  are  shown  in  Fig.  23. 

Footing  i  is  placed  on  the  surface  of  the  layer.  Footing  2  is 
embedded  to  a  depth  of  0.4  x  vq  and  has  a  free  periphery.  Footing  3  is 
embedded  to  the  same  depth  as  Footing  2  but  its  periphery  is  in  contact 
with  the  surrounding  soil.  No  slippage  occurs  where  the  footings  and  the 
layer  are  in  contact. 

The  static  stiffness  coefficients  of  Footing  1,  2  and  3  are 

3 

kg  =  5.34,  7.04  and  11.9  x  G  rQ,  respectively,  where  G  is  the  shear 
modulus  and  rQ  is  the  footing  radius.  This  means  that  the  static 
compliance  is  strongly  reduced  by  the  embedment  of  the  footing  and  by 
side  friction. 

The  above  value  for  the  static  stiffness  of  Footing  3  agrees  well 
with  Kaldjian’s  [21]  results  for  the  torsional  stiffness  of  embedded 
circular  footings. 

The  reduction  of  the  static  compliance  causes  the  resonance 
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frequency  to  increase  and,  as  the  amount  of  geometric  or  radiation 
damping  grows  with  frequency,  this  results  in  an  increased  damping  of 
the  footing  motion.  The  peak  amplitude  is  therefore  even  further  reduced 
by  embedment  of  the  footing  than  is  the  static  compliance. 

A  comparison  of  the  response  curves  for  the  inertia  ratios  Bg=2 
and  Bg=4  shows  that  a  decrease  of  the  inertia  ratio  has  a  similar  effect 
on  the  footing  response  as  have  embedment  and  side  friction. 

The  displacement  functions  from  which  the  response  curves  in  Fig. 

23  were  ca'culated  by  Eqs.  (11.3),  (11.4)  and  (11.6)  are  presented  in 
Fig.  24. 

The  Fj-curves  are  changed  only  insignificantly  by  embedment  and 
side  friction.  The  F^-curve  for  Footing  3  is  somewhat  higher  than  those 
for  Footing  1  and  2.  This  means  that  the  radiation  damping  is  increased 
by  embedment  with  side  friction. 


11.5  Shear  Modulus  Increasing  with  Depth 

In  the  previous  example  it  was  assumed  that  the  medium  supporting 
the  footing  was  homogeneous  with  respect  to  the  mass  density  and  the 
shear  modulus.  However,  the  elastic  moduli  of  soils  usually  vary  with 
depth  even  though  the  mass  density  and  other  properties  of  the  soil  may 
be  fairly  constant  within  the  zone  of  interest. 

The  shear  moduli  of  sands  and  gravels  increase  approximately  as 
the  square  root  of  the  effective  confining  pressure  [43]  and  therefore 
with  the  depth  below  the  surface.  In  the  following  examples,  the  relation 
between  the  shear  modulus  and  the  depth  is  assumed,  accordingly,  tc  be 

G(z)  =  G  /z/r  (11.10) 

o  o 


in  which  z  is  the  depth  below  the  surface,  rQ  is  the  footing  radius  and 

5  is  the  shear  modulus  at  z  =  r  . 
o  o 

The  confining  pressure  in  the  soil  below  the  footing  is  decreased 
by  the  load  of  the  footing.  To  accouut  for  this  fact  in  a  very  simple 
manner  it  is  assumed  that  the  modulus  under  the  footing,  r<ro,  is 

G(z)  =  Gq  */p  for  yz<q  (11.11) 
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in  which  Y  is  the  unit  weight  of  the  soil,  q  is  the  average  contact 
pressure  under  the  footing,  and  p  is  the  dimensionless  load  ratio  de¬ 
fined  by 


«  _SL 


(11.12) 


Fig.  25  shows  the  assumed  variation  of  the  shear  modulus  with  depth  for 

r<r  at  p=0.5,  1  and  2  and  for  r>r  . 
o  o 

Displacement  functions  for  the  torsional  motion  of  three  footings 
without  inertia  are  shown  in  Figs.  27,  28  and  29.  One  footing  sits  on 
the  surface  of  the  ground,  the  other  two  are  embedded  in  the  soil  to  a 
depth  of  0.4  x  r^,  one  with  and  one  without  soil  contact  at  their  sides. 
Each  footing  is  supported  by  a  thick  soil  layer  over  a  rigid  base, 

H/rQ  =  8.  The  shear  modulus  in  the  layer  increases  with  depth  according 
to  Fig.  25  and  its  value  under  the  footing  is  determined  by  the  parameter 
p.  Displacement  functions  for  p  =  0.5,  1.0  and  2.0  are  presented  and 
for  comparison  the  displacement  functions  for  a  homogeneous  layer  are 
also  shown  in  Figs.  27,  28  and  29. 

The  displacement  functions  for  the  layer,  in  which  the  shear 
modulus  increases  with  depth,  differ  somewhat  from  the  displacement 
functions  for  the  homogeneous  layer.  The  difference  is  pronounced  in  the 
case  where  the  footing  sits  on  the  surface  of  the  layer,  but  is  less 
significant  in  the  more  practical  cases  of  embedment  either  with  or 
without  side  friction. 

The  parameter  p  has  hardly  any  influence  on  the  displacement 
functions,  but  it  affects  the  static  spring  constant  kg,  values  of  which 
are  also  presented  in  Figs.  27,  28  and  29  for  the  various  cases. 


11.6  Displacements  in  the  Far  Field 


Displacement  amplitudes  versus  depth  at  different  distances  from 
the  axis  of  symmetry  and  along  the  surface  are  shown  for  three  particular 
cases  in  Figs.  30-33.  In  Case  A  the  footing  sits  on  the  surface  of  a 
homogeneous  layer;  in  Case  B  the  footing  is  embedded  in  the  same  homo¬ 
geneous  layer;  in  Case  C  the  footing  is  embedded  in  a  layer  in  which 
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the  shear  modulus  increases  with  depth  as  shown  in  Fig.  26  (p=1.18). 

In  each  case  the  thickness  of  the  layer  is  equal  to  18  footing  radii, 
the  material  damping  of  the  xayers  is  8  =  3%  of  critical  damping,  and 
the  footings  are  exc  "e.  at  the  frequency  to  =  60tt  so  that  the  amplitudes 
of  the  footing  rotation  are  one  radian. 

Figs.  30  and  31  show  that  the  displacement  amplitudes  close  to  the 
axis  (r>-4ro)  are,  as  expected,  much  larger  at  the  surface  than  at  some 
depth.  But  the  displacement  amplitudes  at  greater  distance  from  the 
axis.  (r=36rQ)  do  not  decrease  with  depth.  This  can  be  explained  by 
t'l  'act  Unfit  u!;e  generalized  Love  waves  generated  by  the  torsional 
vit.i  ions  of  tV «  feltings  are  not  surface  waves,  as  are  the  Rayleigh 
wave.s,  whic  i  h:  amplit  ides  close  to  the  surface  and  attenuate 

with  d  pth,  bn:  •.«.,i.wJ.st  of  S-waves  which  are  reflected  at  the  free  surface 
and  the  rigid  base. 

A  comparison  of  Figs.  30  and  31  shows  that  for  Case  B  the  displace¬ 
ment  ampli .udes  in  the  far  field  are  almost  exactly  twice  as  large  as 
for  case  A.  Tfco  reason  for  this  is  that  the  torque  applied  through  the 
footing  to  the  layer  is  2.04  times  larger  in  Case  P  than  in  Case  A. 

Since  the  wave  length  of  S-waves,  2g,  is  large  compared  to  the  footing 
radius  at  the  considered  frequency,  (&s=8. 6rQ) ,  Saint-Venant's  principle 
applies.  A  different  distribution  of  the  shear  stresses  around  the 
tooting  has  a  negligible  effect  on  the  displacements  at  some  distance 
from  the  footing  if  the  applied  total  torque  is  the  same.  For  high 
frequencies  at  which  the  length  of  S-waves  is  not  long  compared  to  the 
footing  dimension,  this  is  not  true. 

In  Case  C,  Fig.  32,  the  displacement  amplitudes  attenuate  with 
depth  even  at  large  distance  from  the  axis  op  symmetry  (i-36ro) ,  because 
the  shear  modulus  increases  with  depth.  Most  of  the  energy  is  propagated 
through  the  softer  part  of  the  layer  near  the  .rurface. 

This  is  also  the  reason  for  the  amplitudes  of  the  surface  displace¬ 
ments,  Fig.  33,  being  much  larger  in  Case  C  than  in  Case  A  or  Case  B. 
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12.  VERTICAL  MOTION  OF  RIGID  CIRCULAR  FOOTINGS 


12.1  Introduction 

The  vertical  motion  of  a  rigid  circular  footing  which  is  supported 
by  a  homogeneous,  linearly  elastic  half  space.  Fig.  1,  and  excited  by  a 
harmonic  vertical  fcrct.  v?as  first  studied  by  Reissner  [39]  in  1936,  who 
assumed,  for  simplicity.  &  uniform  stress  distribution  under  the  footing. 
Bycroft  [10]  investigated  the  same  problem.  He  assumed  that  the  dynamic 
stress  distribution  beneath  the  rigid  footing  is  similar  to  that  for  the 
static  case  and  obtained  an  approximate  solution  valid  for  low  frequencies. 
Lysmer  [24]  and  Awojobi  and  Grootenhuis  [7]  derived  solutions  for  the 
same  problem  by  accounting  properly  for  the  frequency  dependence  of  the 
stress  distribution  beneath  the  footing. 

The  application  of  the  half  space  theory  to  the  analysis  of  footing- 
soil  systems  implies  the  assumption  that  the  subsoil  under  the  footing 
is  fairly  homogeneous  down  to  a  depth  of  several  footing  diameters. 

Often,  this  assumption  cannot  be  justified,  especially  if  large  footings 
are  to  be  analyzed.  Considerable  deviations  from  the  half  space  theory 
occur  in  cases  in  which  rock  or  hard  layers  are  encountered  at  shallow 
depth  below  the  footing  [29]. 

Warburton  [46]  studied  the  vertical  vibration  of  a  rigid  circular 
footing  supported  by  a  homogeneous  elastic  layer  which  extends  infinitely 
in  the  horizontal  direction  and  rests  on  a  rigid  base.  Fig.  2.  He  assumed 
a  frictionless  interface  between  the  layer  and  the  rigid  base  and  made 
an  approximation  regarding  the  stress  distribution  beneath  the  footing. 
Paul  [38]  studied  the  same  problem  without  assuming  a  priori  a  stress 
distribution  beneath  the  footing.  But  he  did  not  present  a  readily 
useable  solution  and  left  the  results  in  implicit  form. 

In  <:ach  of  the  analytical  half  space  and  layer  solutions  referenced 
above,  it  was  assumed  that  the  interface  between  the  footing  and  the 
supporting  medium  is  frictionless. 
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Lysmer  and  Kuhlemeyer  [25]  used  successfully  a  finite  element 
method  with  an  energy  absorbing  boundary  (see  Chapter  1)  to  study  the 
effect  of  embedment  on  the  vertical  motion  of  a  footing  which  is  supported 
by  a  homogeneous  elastic  half  space.  Kuhlmeyer  [22]  used  the  same  method 
to  study  the  vertical  motion  of  a  footing  supported  by  a  stratified  half 
space  and  found  that  the  method  gave  inconclusive  results  when  the  deeper 
layers  were  stiffer  than  the  surface  layer,  which  is  usually  the  case 
in  foundation  vibration  problems. 

In  the  following  examples  the  vertical  motion  of  rigid  circular 
footings  on  the  surface  of  a  homogeneous  layer  and  of  those  embedded  in  a 
layer  in  which  the  elastic  moduli  increase  with  depth  are  investigated 
by  the  method  presented  in  this  dissertation. 

Before  the  examples  can  be  discussed  a  few  terms  need  to  be  ex¬ 
plained.  The  harmonic  vertical  motion  of  a  rigid,  massless  footing  is 
described  by 

P 

w  =  —  F  exp(icjt)  (12.1) 

z 

in  whi^h  w  is  the  vertical  displacement  of  the  footing,  Pq  is  the 
amplitude  of  the  vertical  force  exciting  the  footing  motion  and  kz  is  the 
static  spring  constant.  The  dimensionless  complex  quantity  F  =  F.^  +  i  F2 
is  the  displacement  function  as  before  and  depends  on  the  circular 
frequency  u)  or  the  dimensionless  frequency  ratio  aQ  defined  by  Eq.  11.2. 

If  the  displacement  function  F  for  a  massless  footing  is  known, 

one  can  calculate  the  response  of  a  footing  with  the  mass  m,  as  shown 

by  Lysmer  and  Richart  [24],  by 

P 

w  =  ~  M  cos(u)t  +  <j>)  (12.2) 


in  which  M  is  the  dynamic  magnification  factor 


M  = 


2  2 

F*  +  F 
1  2 


(1-u^m  F,/k  )^  +  (u^m  F0/k  )^ 

I  2  x.  2 


(12.3) 


and  <j>  is  the  phase  shift 


<)i  =  tan 


Fj-c^m  (F2  +  F2)/kz 


(12. A) 


A  dimensionless  mass  ratio  may  be  defined  by 

G  r 

„  m  o 

b2  ■ — 3 '  — 

p  r  2 


(12.5) 


in  which  p  and  G  are  the  mass  density  and  the  shear  modulus  of  the  soil, 
respectively,  and  rQ  is  the  footing  radius. 


12.2  Footing  on  Homogeneous  Layer 

The  vertical  motion  of  a  rigid  circular  footing  supported  by  a 
homogeneous  layer.  Fig.  2,  is  studied.  The  footing  has  the  radius  rQ  and 
is  smooth  at  its  interface  with  the  layer,  permitting  free  hori2ontal 
displacements  of  the  layer  surface  beneath  the  footing.  The  layer  has 
the  thickness  H  and  is  welded  to  the  rigid  base.  Poisson's  ratio  of  the 
layer  material  is  V  =  0.25  and  the  complex  shear  modulus  G=(l  +  i  2  3) 
units,  where  8  is  the  fraction  of  critical  damping. 

Displacement  functions  for  the  ratio  of  layer  thickness  to  footing 

radius  H/rQ  =  10  and  6  are  presented  in  Figs.  34  and  35  for  different 

values  of  8.  The  curves  for  H/r  =  10  and  H/r  =  6  look  similar  if  the 

o  o 

abscissas  in  Figs.  34  and  35  are  scaled  by  the  layer  thicknesses.  In 
the  undamped  case,  6=0%,  the  displacement  functions  show  singularities, 
but  in  the  damped  case,  6=5%  or  10%,  they  are  "smooth." 

To  explain  the  singularities  in  the  undamped  case  it  is  best  to 
look  at  the  spectral  lines  of  the  wave  numbers  which  characterizes  the 
wave  modes  used  in  the  displacement  expansion  according  to  Eq.  (7.9). 

The  first  five  spectral  lines  are  presented  in  Fig.  3b  which  shows 
a  three-dimensional  space.  The  axes  in  Fig.  36  are  labeled  £=Fe(k)*H, 
0=Im(k)*H  and  8=uH/Vg  where  II  is  the  thickness  of  the  layer.  The  long 
dashed  lines  lie  in  the  plane  and  represent  real  wave  numbers.  The 

short  dashed  lines  lie  in  the  q-Q-  plane  and  represent  imaginary  wave 
numbers.  The  full  lines  are  curves  in  the  5-tj-Q-  space  and  occur  in 
pairs  which  are  symmetric  about  the  ri-ft-  plane;  each  pair  represents  a 


pair  of  complex  wave  numbers,  k  =  kj-ik2  and  -k  =  -kj-ik,,. 

The  spectral  lines  intersect  the  planes  0=0,  5=0  and  ri=0  at  right 

angles.  This  fellows  from  the  fact  that  the  secular  equation,  Eq.  6.27, 

2  2 

is  an  analytical  function  of  u>  and  k  in  the  undamped  case.  A  proof  can 
be  performed  in  a  similar  manner  to  that  given  by  McNiven  et  al.  [31]  for 
the  spectral  lines  of  waves  in  axisymmetric  rods  and  will  therefore  not 
be  presented. 

The  points  at  which  the  spectral  lines  intersect  the  0-axis  define 
natural  frequencies  of  the  layer  corresponding  to  vertically  propagating 
waves  which  are  reflected  at  the  free  surface  and  the  rigid  base  resulting 
in  a  standing  wave.  At  these  frequencies  the  dsiplacements  in  the  layer 
vary  only  with  the  depth  z  and  any  plane  parallel  to  the  surface  remains 
plane  during  motion  because  k=0. 

The  natural  frequencies  labelled  a,  c  and  d  in  Fig.  36  are 
associated  with  pure  S-waves  producing  only  horizontal  displacements. 

They  can  be  computed  by 

wn  =  *  Vs/H-  n=l,2,  ...  (12.6) 


in  which  is  the  shear  wave  velocity.  Because  the  waves  corresponding 
to  the  points  a,  c  and  d  have  no  displacement  components  in  the  vertical 
direction,  they  do  not  affect  the  vertical  motion  of  the  footing. 

The  natural  frequencies  labelled  b  and  e  in  Fig.  36  are  associated 
with  pure  P-waves  producing  only  vertical  displacements  in  the  layer. 

The  layer  surface  can  move  up  and  down  at  these  frequencies  without  being 
ex-.ited  by  external  forces,  and  since  it  remains  plane  the  boundary 
conditions  under  the  massless  footing  are  identical  to  those  of  the  free 
surface.  Therefore,  these  frequencies  are  natural  frequencies  of  the 
footing-layer  system  if  the  footing  is  massless.  Consequently,  the  dis¬ 
placement  functions  for  3=0  in  Figs.  34  and  35  have  singularities  at  the 
corresponding  frequency  ratios.  The  natural  frequencies  in  question  are 


"a  *  Hf1  *  VH 


n— 1 , 2 ,  . . ■ 


(12.7) 


in  which  V  is  the  velocity  of  P-waves. 

The  first  two  natural  frequencies  for  the  vertical  motion  of  the 


massless  footing  in  the  case  H/r  =  10  occur  at  a  =0.272  and  a  =0.816. 

o  oo 


Fig.  34  shows  clearly  the  singularity  at  a  =0.272,  but  the  singularity 
at  ao=0.8l6  is  hardly  noticeable. 

In  Fig.  34  one  observes  another  irregularity  at  a^=Q. 746  which 
corresponds  to  the  frequency  at  g  and  g'  in  Fig.  36.  At  this  frequency 
two  wave  modes  exist  in  the  layer  which  have  equal  and  opposite  phase 
velocities.  The  group  velocity,  which  is 


(12.8) 


if  k  is  real  [16],  vanishes  at  g  and  g'.  The  two  wave  modes  have, 
therefore,  zero  group  velocities  at  this  frequency  and  do  not  transmit 
energy,  because  the  rate  of  energy  transmitted  by  a  wave  is  the  product 
of  the  group  velocity  and  the  energy  density  in  the  wave  [17]. 

If  the  amplitudes  of  the  two  wave  modes  are  equal,  the  modes 
combine  to  a  standing  wave  and  the  horizontal  and  vertical  displacements 
in  the  axisymmetric  case  vary  along  the  free  surface  as 


<5x  «  (kr)  exp(i/idt)  (12.9a) 

6z  «  JQ  (kr)  exp(iuit)  (12.9b) 


in  which  J  and  J,  are  the  Bessel  functions  of  order  zero  and  one, 
o  1 

respectively.  Eq.  (12.9)  follows  from  Eqs.  (6.32,  7.1,  7.3)  and  from 
the  properties  of  the  Hankel  functions  [2]. 


H(2,(kr)  -  H(2)(-kr)  =  2J  (kr)  (12.10a) 

o  o  o 

H*2)(kr)  +H^2)(-kr)  =  2^  (kr)  (12.10b) 


Because  cne  wave  number  corresponding  to  point  g  in  Fig.  36  is 

relatively  small,  kr  =0.208  for  H/r  =  10,  the  free  surface  in  the  area 

o  o 

r<ry  remains  approximately  pla  ie  during  vibration  in  these  two  modes. 

The  difference  in  the  vertical  displacements  at  r=rQ  and  r=0  is  only  one 
percent.  This  means  that  the  displacements  of  the  free  surface  are 
almost  compatible  with  those  under  the  rigid  massless  footing.  Therefore 
it  takes  only  a  small  force  to  completely  satisfy  the  boundary  conditions 
under  the  footing  and  to  produce  a  large  footing  displacement  at  this 
frequency.  Consequently,  the  F^-curve  for  8=0 7.  in  Fig.  34  shows  a  peak 
at  aQ=0.746.  The  peak  is  finite  in  height;  but  it  is  difficult  to 
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|  determine  the  height  numerically,  because  the  peak  is  very  narrow. 

k 

f  The  same  phenomenon  occurs  at  the  frequency  corresponding  to  the 

|  points  f  and  f*  in  Fig.  36.  But  since  this  frequency  is  very  close  to 

*  the  lowest  natural  frequency  of  the  footing- layer  system,  no  separate 

I  peak  shows  up  in  Fig.  34. 

[  A  similar  phenomenon  to  that  described  above  was  observed  experi- 

j  mentally  in  steel  rods  by  Oliver  [36]  and  studied  theoretically  by 

;  McNiven  [30).  In  the  rods,  very  narrow  and  high  response  peaks  occurred 

i 

also  at  those  frequencies  at  which  pairs  of  complex  spectral  lines  for 
wive  numbers  became  purely  real  as  happens  at  the  points  f  and  f*  or  g 
and  g'  in  Fig.  36. 

The  effect  of  material  damping  on  the  response  of  footings  is  to 
smooth  the  displacement  functions  and  eliminate  natural  frequencies. 

,  lhe  displacement  functions  in  Fig.  34  and  35  for  8=5%  and  10%  of 

1 

;  critical  damping  therefore  show  no  singularities. 

t 

!  The  spectral  lines  for  the  wave  numbers  in  the  damped  case  are  not 

j  much  different  from  those  shown  in  Fig.  36  for  the  undamped  case.  But 

l  they  do  not  intersect  the  plane,  as  all  wave  numbers  have  negative, 

(  non-zero  imaginary  parts  in  the  damped  case,  and  they  are  smooth  with 

! 

!  continuous  first  derivatives  at  the  frequencies  at  which  the  spectral 

lines  form  right  angles  in  the  undamped  case. 

It  is  more  realistic  to  include  some  material  damping  in  the 
analysis  than  to  assume  that  the  soil  behaves  purely  elastically. 
Furthermore,  the  incorporation  of  damping  has  the  advantage  that  singular¬ 
ities  are  avoided. 

Response  spectra  for  footings  with  different  mass  ratios  are 
presented  in  Figs.  37  and  38  for  H/r^  =  6  and  H/r^  =  10  at  8-5%.  They 
show  peaks  at  frequencies  close  to  the  first  natural  frequencies  of  the 
corresponding  undamped  cases.  The  peaks  are  high  for  H/rQ  =  6,  but  not 
very  significant  for  H/rQ  =  10. 

Except  for  these  peaks,  the  response  spectra  for  H/r^  =  10  are 

quite  similar  to  those  for  the  vertical  motion  of  a  rigid  circular 

footing  on  a  homogeneous  elastic  half  space,  which  are  shown  in  Fig.  39 

and  are  based  on  Lysmer's  solution  [24] .  However,  the  response  spectra 

for  H/r  =  6  differ  markedly  from  those  for  the  half  space, 
o 
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The  values  of  the  displacement  functions  in  Figs.  34  and  35  were 

computed  by  using  rectangular  meshes,  which  were  similar  to  thqt’  shown 

in  Fig.  21  and  consisted  of  22  x  5  finite  elements  and  22  sublayers. 

The  maximum  dimension  of  any  element  and  the  maximum  thickness  of  any 

layer  was  less  than  1/10  of  the  length  of  shear  waves  at  the  highest  i 

frequency  considered.  v  N 

To  verify  the  validity  of  the  theory  and  to  check  the  computer 

program,  analyses  for  H/rQ  =  6  were  repeated  at  several  frequencies  with 

a  different  mesh,  in  which  a  much  larger  region,  r<6rQ  as  compared  to 

r<rc,  was  represented  by  finite  elements.  The  results  obtained  using 

the  two  different  meshes  were  virtually  identical.  The  relative 

differences  between  the  computed  values  of  the  displacement  functions 

( 

were  less  than  one  percent. 


12.3  Footing  Embedded  in  Inhomogeneous  Layer 


Fig.  40  shows  an  embedded  rigid  circular  footing  which  is  forced 
to  vibiate  in  the  vertical  mode.  The  shear  modulus  in  the  supporting 
layer  increases  with  depth  as  defined  in  Fig.  26.  Poisson’s  ratio  is 
v=l/3  throughout  the  layer,  and  3%  of  critical  damping  for  sh'-ar  and 
dilatational  deformations  in  the  layer  is  assumed.  The  damping  is 
introduced  by  using  a  complex  shear  mor/'ius  (Fig.  26  shows  only  the  real 
part  of  the  shear  modulus).  Since  Poisson’s  ratio  is  real.  Lame's 
constant  X=2Gv/ (1-2V)  is  also  complex.  The  layer  has  c.  thickness  of 
18  footing  radii  and  rests  on  a  rough  rigid  base. 

The  displacement  function  for  a  massless  footing  and  the  response 
spectrum  f.ir  a  footing  with  mass  are  presented  in  Fig.  40.  The  figure 
indicates  a  small  peak  at  a  frequency  close  to  that  at  which  the  layer, 
in  the  absence  of  damping,  would  have  its  first  natural  frequency  for 
vertical  displacements. 

A  comparison  with  the  response  spectra  in  Figs.  37  and  38  suggests 
that  the  response  peak  in  Fig.  40  is  small  because  the  layer  is  thick 
with  respect  to  t.he  footing  radius,  H/ rQ  =  18.  Due  to  damping,  the 
displacement  function  shows  no  further  peaks  at  the  frequencies  at 
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which  the  layer,  in  the  absence  of  damping,  would  have  higher  natural 
frequencies. 

Figs.  41  and  42  show  displacement  amplitudes  along  the  surface  of 
the  layer  and  displacement  amplitudes  versus  depth  at  various  distances 
from  the  axis  of  the  footing.  The  displacements  are  caused  by  vertical 
vibration  of  the  footing  with  unit  displacement  amplitude.  The  horizontal 
and  vertical  components  of  the  displacements  along  the  surface  have 
about  equal  amplitudes  except  close  to  the  foocing  where  the  vertical 
component  dominates.  Fig.  42  shows  that  the  displacement  amplitudes 
decrease  with  depths  at  close  and  far  distances  from  the  axis  of  the 
footing.  This  means  that  most  of  the  energy  is  propagated  in  the  upF-r 
part  of  the  layer,  where  the  material  is  softer  chan  at  greater  dep'.h, 
and  also  that  the  rigid  base  assumed  in  the  analyses  has  little  or  no 
influence  on  the  results. 


13.  SUMMARY  AND  CONCLUSIONS 


A  numerical  method  has  been  presented  for  the  analysis  of  steady- 
state  wave  propagation  problems  in  linearly  elastic  or  viscoelastic 
media  of  infinite  extent.  The  geometries  considered  are  either  plane  or 
axisymmetric  and  consist  of  a  finite,  irregular  region  which  is  joined 
to  semi-infinite  layered  regions  of  the  type  shown  in  Figs.  4  and  5.  All 
external  loads  are  applied  within  the  irregular  region  and  vary  harmonic¬ 
ally  with  time.  Motions  both  in  the  plane  and  perpendicular  to  the 
plane  of  a  cross-section  through  an  axisymmetric  or  plane  geometry  have 
been  treated. 

The  irregular  region  is  discretized  by  the  finite  element  method. 
Compatible  finite  elements  with  quadrilateral  cross-sections  are  employ* 
and  the  displacements  at  the  nodal  joints  of  the  elements  are  introduced 
as  the  degrees  of  freedom. 

The  semi- infinite  iayered  region  is  discretized  by  subdividing  the 
layers  into  thin  sublayers  and  by  assuming  that  within  each  sublayer  the 
displacements  vary  linearly  in  the  direction  normal  to  the  layers.  In 
the  direction  parallel  to  the  layers,  the  displacements  are  required  to 
satisfy  the  pertinent  ordir.  iry  differential  equations,  which  are  obtained 
by  separation  of  the  variables  from  the  partial  differential  equations 
that  govern  the  motion  of  the  continuum. 

Compatibility  conditions  for  displacements  and  equilibrium  must  be 
satisfied  along  the  sublayer  interfaces,  the  stress  free  surface  and 
the  fixed  base.  These  conditions  constitute  an  algebraic  eigenvalue 
problem,  since  no  external  forces  are  applied  within  the  layered  region. 
The  solutions  to  the  eigenvalue  problem  describe  the  propagating  and 
standing  waves  which  can  exist  in  the  discretized  layered  region.  These 
waves  serve  as  shape  functions  for  the  displacement  expansion.  A  dynamic 
stiffness  matrix  is  then  developed  which  uniquely  relates  the  nodal 
forces  to  the  simultaneous  nodal  displacements  and  properly  represents 
the  elastic  ana  viscous  dynamic  response  of  the  semi-infinite  layered 
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region. 

The  equations  of  motion  for  the  discrete  model  are  derived  from 
the  principle  of  virtual  work  which  has  been  formulated  to  facilitate 
the  use  of  complex  variables  in  the  analysis  of  harmonic  motion. 

An  important  step  in  the  procedure  is  the  determination  of  the 

wave  modes  in  the  layered  region  by  the  solution  of  the  eigenvalue 

problem.  For  motions  perpendicular  to  the  cross-sectional  plane  the 

problem  ([A]k  +  f C] )  {v}  -  {Ol  must  be  solved  in  which  [A]  and  [C]  are 

complex,  symmetric,  tridiagonal  matrices.  For  motions  in  the  cross- 

2 

sectional  plane  the  problem  is  ([A]k  +  i[B]k  +  [C])  {v}  =  {0}  in  which 
f A j  and  [C]  are  complex,  symmetric  matrices  with  bandwidth  5  and  [B] 
is  a  complex,  skew-symmetric  matrix  with  bandwidth  7.  Efficient  FORTRAN 
IV  subroutines  for  solving  these  problems  have  been  developed  and 
presented. 

The  method  is  consistent  in  that  the  same  displacement  expansions 
are  employed  for  the  derivation  of  the  elastic,  damping  and  inertia 
forces  at  any  nodal  joint.  The  displacement  expansions  both  in  the 
finite  element  regions  and  in  the  layered  regions  are  compatible  and 
thus  provide  continuous  displacements  throughout  the  entire  domain.  The 
analysis  yields  the  displacements  in  the  irregular  region  as  well  as  in 
the  semi- infinite  layered  regions. 


Accuracy  Studies 

The  accuracy  of  the  method  was  investigated  for  the  example  in 
which  a  line  load  acts  on  the  surface  of  a  homogeneous  layer  and 
generates  plane  generalized  Love  or  Rayleigh  waves.  In  the  Love  wave 
case,  the  numerical  results  for  the  discretized  layer  were  compared  with 
the  exact  solution  for  the  continuous  layer.  Good  agreement  was  obtained 
with  a  relatively  coarse  mesh  and  excellent  agreement,  was  achieved  with 
a  fine  mesh.  In  the  Rayleigh  wave  case,  the  accuracy  of  the  numerical 
solutions  nad  to  be  -tudied  indirectly  by  comparing  solutions  obtained 
with  different  meshes,  because  no  exact  continuum  solution  is  available. 
For  this  purpose  the  boundary  between  the  irregular  and  the  layered 


region  was  shifted,  replacing  the  finite  element  displacement  field  by 
the  expansion  into  wave  modes  and  vice  versa.  The  size  of  the  finite 
elements  and  the  thickness  of  the  sublayers  were  also  varied.  The 
relative  differences  between  numerical  solutions  obtained  with  different 
meshes  were  insignificant.  In  addition,  exact  wave  numbers  for  waves 
in  the  continuous  layer  were  computed  and  compared  with  the  wave  numbers 
for  the  corresponding  waves  in  the  discretized  layer.  The  agreement  was 
highly  satisfactory  and  improved  as  the  discretization  was  refined. 
Similar  studies  of  numerical  solutions  to  axisymmetric  problems  were 
performed  but  have  not  been  presented  here.  In  these  cases  the  agreement 
of  solutions  obtained  with  different  meshes  was  equally  good. 

In  another  study,  a  numerical  solution  for  the  torsional  vibration 
of  a  circular  footing  on  a  homogeneous  layer  was  computed  for  comparison 
with  published  analytical  solutions  and  test  results.  The  numerical 
method  was  again  found  to  yield  excellent  results. 


Torsional  Motion  of  Circular  Footing 


The  method  presented  herein  was  primarily  developed  for  the  analysis 
of  steady-stat  vibrations  of  machine  foundations.  Thus  emphasis  was 
placed  on  applications  to  this  type  of  problem. 

The  torsional  response  of  circular  footings  to  vibratory  torques 
about  the  vertical  axes  of  the  footings  was  studied.  The  effects  of  the 
thickness  of  a  supporting  layer,  embedment  of  the  footing  and  increasing 
shear  modulus  with  depth  were  investigated.  The  response  of  a  footing  on 
a  homogeneous  layer  which  has  a  thickness  of  4  footing  radii  .'r  more  is., 
for  practical  purposes,  identical  to  that  of  a  footing  on  a  homogeneous 
half  space.  Footings  on  layers  with  a  thickness  of  one  footing  radius 
or  less  respond  quite  differently  and  do  not  benefit  from  radiation 
damping  at  low  frequencies. 

Embedment  may  reduce  the  vibration  amplitudes  drastically, 
especially  when  the  footing  sides  are  bonded  to  the  adjacent  soil. 

Several  cases  were  studied  in  which  the  shear  modulus  increased  as  the 
square  root  of  the  depth  below  the  ground  surface.  A  careful  determination 
of  the  static  stiffness  is,  as  usual,  of  primary  importance.  The  dynamic 
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effects  are  similar  to  those  in  the  case  of  a  homogeneous  half  space. 
Therefore,  the  displacement  functions  of  the  half  space  solution  may  be 
utilized  in  the  analysis  as  a  first  approximation  provided  that  the  shear 
modulus  at  a  depth  of  about  one  half  of  the  footing  radius  is  chosen  as 
the  representative  shear  modulus.  A  study  of  displacements  in  the  far 
field  has  indicated  that,  in  the  typical  frequency  range,  the  amplitudes 
of  surface  displacements  at  some  distance  from  the  footing  depend  pre¬ 
dominantly  on  the  torque  transmitted  into  the  ground  and  on  the  layering 
of  the  soil. 


Vertical  Motion  of  Circular  Footings 
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The  vertical  vibration  of  a  circular  footing  supported  by  a  homo¬ 
geneous  layer  over  a  rigid  base  was  investigated  for  two  ratios  of  layer 
thickness  to  footing  radius,  H/rQ  =  6  and  H/rQ  =  10.  Displacement 
functions  of  massless  footings  were  computed  for  cases  with  and  without 
damping.  In  the  undamped  case  two  different  resonance  phenomena  were 
observed  and  discussed.  Response  spectra  for  footings  with  typical  mass 
ratios  show  only  cne  pronounced  peek  when  a  realistic  amount  of  damping 
is  assumed.  This  peak  occurs  at  the  fundamental  natural  frequency  for 
vertically  travelling  P-waves  in  the  soil  layer.  For  H/r^  =  6  or  smaller 
this  peak  is  high,  but  for  H/rQ  =  10  or  larger  this  peak  is  insignificant. 
The  response  of  footings  on  a  thick  layer,  H/rQ  >  10,  is  similar  to  that 
of  a  footing  on  a  homogeneous  half  space.  The  vertical  response  of  a 
circular  footing  embedded  in  a  layer,  in  which  the  moduli  increase  with 
depth,  was  also  analyzed. 


Screening  Effect  of  Trenches 

The  versatility  of  the  method  was  demonstrated  by  a  further 
example.  The  screening  effect  of  trenches  on  horizontally  polarized 
S-waves  was  investigated  by  varying  the  depth  of  the  trench  and  the 
frequency  of  the  wave  motion.  It  was  shown  that  trerches  can  be  quite 
effective  as  a  barrier  against  waves  of  some  frequencies  but  may  be 
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ineffective  for  waves  of  other  frequencies. 


Possible  Additional  Applications 

It  is  believed  that  the  method  may  also  be  of  value  in  the  study  of 
soil-structure  interaction  during  earthquakes.  This  would  require  the 
introduction  of  a  few  features  in  addition  to  those  included  in  the 
analyses  presented  herein.  In  this  type  of  problem,  the  rigid  base  of 
the  model  may  be  subjected  to  motion  in  the  form  of  a  time-displacement 
history  of  the  bedrock  underlying  the  soil.  The  induced  transient  motion 
may  be  analyzed  with  the  assistance  of  a  discretized  Fourier  transforma¬ 
tion. 

The  method  might  also  find  application  to  wave  propagation  problems 
in  other  fields.  In  seismology,  fault  mechanisms  of  earthquakes  and  the 
effect  of  valleys  and  mountain  ranges  on  seismic  waves  may  be  studied. 

The  interpreation  of  results  from  non-destructive  vibratory  pavement 
testing  may  also  be  improved  with  the  aid  of  this  method  of  analysis. 
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FIG.  6-DISCRETE  AXISYMMETRIC  MODEL 


FIG.  7-  HYSTERETIC  STRESS -STRAIN  RELATIONSHIP 
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FIG.  8- REPRESENTATION  OF  STRESS  AND  STRAIN 
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FIG.  16**  LINE  LOAD  ON  HOMOGENEOUS  LAYER 
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FIG.  25-  SHEAR  MOOULUS  INCREASING  WITH  DEPTH 
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FIG. 26 -SHEAR  MODULUS  VS.  DEPTH  FOR  INHOMOGENEOUS  LAYER 
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FIG  32- DISPLACEMENT  AMPLITUDE  VS.  DEPTH  FOR  TORSIONAL  MOTION 
OF  EMBEDDED  CIRCULAR  FOOTING  ON  INHOMOGENEOUS  LAYER 
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iG.  39-RESPONSE  SPECTRA  FOR  VERTICAL  MOTION  OF  CIRCULAR  FOOTINGS  ON  HALF  SPACE 


FIG.  40- DISPLACEMENT  FUNCTION  AND  MAGNIFICATION  FACTOR  FOR 

VERTICAL  MOTION  OF  CIRCULAR  FOOTING  ON  INHOMOGENEOUS  LAYER 


FIG. 4 1  -DISPLACEMENT  AMPLITUDES  ALONG  SURFACE  FOR  VERTICAL  MOTION  OF 
CIRCULAR  FOOTING  ON  INHOMOGENEOUS  LAYER  AT^=60tt 
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APPENDIX  1 


SOLUTION  METHOD  FOR  ([A]k  +  [C])  {v}  =  {0} 


The  n  x  n  matrices  [A]  and  [C]  in 
( [A]k2  +  [C])  {v}  =  {0} 


(Al.l) 


are  tridiagonal,  complex  symmetric  and  [A]  is  well  conditioned.  The 

diagonal  elements  of  [A]  and  [Cj  are  called  and  c^ ,  j=l,  ...»  n,  and 

the  off-diagonal  elements  are  called  and  d  ^ ,  j=l,  . ..,  n-1,  respectively. 

The  eigenvalues  and  eigenvectors  of  Eq.  (Al.l)  are  generally  complex. 

A  determinant  search  technique  employing  Newton's  iteration 

2 

method  is  used  to  compute  the  n  eigenvalues  X  =  k  ,  s=l,  ...»  n,  of 

s  s 

Eq.  (Al.l). 

The  characteristic  polynomial  to  Eq.  (Al.l),  deflated  by  the 
divisors  which  are  already  known  is 


fr  =  I  [A]  tr+  £C]|*(V(tr-  Xs))  1 


(A1.2) 


in  which  t  is  the  rth  approximation  to  the  mth  root,  X^,  and  .\g, 
s=l,  ...,  m-1,  are  the  roots  already  calculated. 

Newton's  iteration  schema  is 


t  ..  =  t  -  f  /f 
r+1  r  r  r 


(A1.3) 


in  which  *  denotes  the  derivative  with  respect  to  t.  The  diagonal 
elements  after  triangular  factorization  of  ([A] tr  +  [C] )  are 


U1  =  Vr  +  C1 

uj+l  =  aj+ltr  +  cj+l  "  (bjCr  +  V  /uj 


j”l )  •  *  •  j  n**l 


(Ai.4) 


f  -  n  «  /  n  (tr  -  x  ) 

j=i  y  s=i 


(A1.5) 
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and 


The  sth  eigenvector  corresponding  to  the  eigenvalue  A  is  computed 

s 

by  two  steps  of  inverse  iteration 

<IAJAS  +  [ClHVg}^  -  {vs)r  ;  r=l,2  (A1.8) 

where  {v  }.  is  chosen  as  a  full  unit  vector, 
s  1 

A  FORTRAN  IV  subroutine,  ALAMDAB,  which  is  based  on  the  method  out¬ 
lined  above  is  listed  in  Appendix  5.  About  10  iterations  are  needed  on 
the  average  to  obtain  12  to  14  decimal  digits  for  each  eigenvalue. 
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AVp^^3,  **' 


are  symmetric  2m  x  2m  matrices  and 


{w}  = 


W, 

w. 

1 

_ 

1 

W2 

kw^ 

(A2.8) 


The  2m  eigenvalues  of  Eq.  (A2.4)  and  their  corresponding  eigen¬ 
vectors  are  in  general  complex.  The  row  and  column  eigenvectors  are 
identical  because  [E]  and  £F]  are  symmetric.  The  eigenvectors  are 
orthogonal  to  each  other  with  respect  to  [F]  and  can  be  normalized  in 
the  sense 


(w)T  [F]  (w)  =0  if  r^s 

S=1  if  r=s  <A2‘9> 

The  generalized  Rayleigh  quotient  iteration,  GRQI,  is  ch  sen  as 
solution  method  because  it  can  take  full  advantage  of  the  narrow  band- 
widths  of  [A],  [B],  and  [C]  and  it  converges  rapidly.  Ostrowski  [37] 
proved  that,  under  certain  conditions,  the  GRQI's  convergence  rate  is 
cubic  in  the  limit.  GRQI  consists  of  inverse  iteration  with  a  shift  by 
the  generalized  Rayleigh  quotient,  p,  at  each  step 


UE]  -  Pj.ilFD  =  IF] 


(A2.10) 


{x}.  [E]  tx). 


{xK  [F]  {x}. 


P  -  - J - JL  -  p  +  - J - J- 

3  {x}*  [F]  {x}j  j*1  {x}*  [F]  (x}j 


(A2.ll) 


where  {xj^  [F]  {x)j  ^  0.  As  j  °°,{x}j  and  will  go  to  (w)  and  k, 
respectively.  The  generalized  Rayleigh  Quotient,  Eq.  (A2.ll),  is  a 
complex  number  and  is  stationary  in  the  neighborhood  of  an  eigenvector 
[37]. 

The  iteration  scheme,  Eqs.  (A2.10)  and  (A2.1I),  is  now  specialized 
to  take  advantage  of  the  structures  of  [E]  and  [F].  The  steps  are 
(i)  Inverse  iteration 

(p  1  [A]  +  +  [C]){xPj=[C]{x1}j_1-pj_1[A]{x2}j_1 

(A2.12) 
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(xi.}j  (txPj  "  txi}j-i)/pj=i 


(A2.13) 


in  which  the  '  indicates  a  vector  which  is  not  normalized; 
Rayleigh  Quotient 


Y?  =  1°]  ^x{^j  ~  ^xpT  lAJ  ^xPj 

Pj  =  pj-l  +  ({xl}j  tCl  {xl}j-l  "  {x2} j  lAl  {x2}j-l)/Yj 


(A2.14) 


(iii)  Normalization  of  vector 

txi}j  =  *xiVpJ 


(A2.15) 


(A2.16) 


(x_}4  =  {x;}./p. 

j  *•  j  i 


(A2.17) 


(iv)  Check 


,P1  -  pJ-ll 


<  tolerance 


then  let  k=p.  and  {w}  - 
J 


otherwise  j=j+l  and  return  to  (i) . 

The  matrix  in  the  lefthand  side  of  Eq.  (A2.12)  has  half  bar- ..width 
4,  including  the  diagonal.  The  triangular  factorization  is  ther  ' .-e  very 
fast.  The  FORTRAN  IV  subroutine  listed  in  Appendix  6  involves  approxi¬ 
mately  40  x  m  complex  multiplications  and  additions  per  iteration.  On 
the  average  7  iterations  were  needed  to  find  about  12  correct  digits 
of  an  eigenvalue  and  its  corresponding  eigenvector  from  an  arbitrary 
initial  vector. 

To  avoid  convergence  towards  already  known  eigenvectors  the 
iteration  vector  is  deflated  by  Gram-Schmidt's  orthogonalization  process 
[49],  In  order  to  keep  the  computational  work  involved  in  the  deflation 
small  an  arbitrary  vector  {p}  of  dimension  2ra  is  selected  which  is 
assumed  to  be  complete  in  the  2m-dimensional  space,  i.e. 


(A2.18) 


{p} 


a  *0 
s 


The  starting  vector  {x}q  for  the  iteration  towards  the  rth  eigenvector 
is  then  obtained  by 

r-1 

(x>  =  {p}  -  y  (w>  a  (A2.19) 

O  L-i  s  s 

s=l 


with 


a  =  {p}T  [F]  (w>  (A2.20) 

s  s 

From  the  structures  of  [A],  [B],  and  [C]  it  follows  that  if 
(k,  {w})  is  a  solution  to  Eq.  (A2.4)  then  (-k,  {w})  is  another  solution 
in  which 


Therefore  only  m  eigensolutions  to  Eq.  (A2.4)  need  to  be  determined  by 
iteration  and  the  deflation  by  Eq.  (A2.19)  can  be  performed  for  a  pair 

of  eigenvectors,  ^^s-l  an(*  ^w^2s  =  ^^2s-l*  simultaneously* 

As  Gram-Schmidt's  orthogonalization  process  is  sensitive  to 
roundoff  errors,  the  eigenvalues  and  eigenvectors  have  to  be  precise. 

No  difficulties  in  finding  all  the  eigenvalues  and  eigenvectors  by  the 
presented  method  have  been  experienced  in  a  few  hundred  analyses  of 
problems  with  dimensions  up  to  2m  =  140  on  the  CDC  6400  computer,  which 
has  a  mantissa  of  48  binary  bits  in  floating  point  arithmetic. 

The  complete  solution  of  a  problem  with  the  dimension  2m  =  80  takes 
about  16  seconds  on  the  CDC  6400.  The  computation  time  increases  with 
the  square  of  the  dimension. 
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APPENDIX  3 


EXACT  SOLUTION  FOR  LINE  LOAD  ON  HOMOGENEOUS  LAYER 


The  physical  system  considered  is  shown  in  Fig.  16.  The  loading 
consists  of  a  harmonic  line  load,  P  =  1  •  exp(iajt)  per  unit  length. 

The  motion  of  the  elastic  layer  is  governed  by  the  wave  equation 


ju  G  /3u  3u\ 

.2  oL.2  .2/ 


(A3.1) 


at  p  3x  3z 

in  which  u  is  the  displacement  in  the  y-direct.ion,  G  is  the  shear 
modulus,  and  p  is  the  mass  density.  Eq.  (A3.1)  has  the  plane  wave 
solution 

u(x,z,t)  =  (A  cos  yz  +  B  sin  yz)  exp(iwt-ikx)  (A3. 2) 

in  which  A  and  B  are  arbitrary  constants,  k  is  the  wave  number,  and 


_  k2 
G  * 


(A3. 3) 


Satisfaction  of  the  boundary  conditions  u(x,H,t)  =  0  and  t  (x^0,z=0,t) 

zy 

=  0  leads  to  the  eigenfunctions 


v  (x,z,t)  =  cos(y  z)  exp(iait-ik  x) 

5  S  S 


in  which 


2s-l  IT 
Ys  =  2  H 


s=l,2. 


(A3. 4) 


(A3. 5) 


,2  _  t/p  r2s-l  tt_.  2 
s  ~  G  ~  2 


(A3. 6) 


An  interpretation  of  the  solutions  defined  by  Eqs.  (A3. 4)  to  (A3. 6) 
and  rules  for  the  choice  of  the  sign  on  kg  have  been  given  previously  in 
Chapter  4. 

The  shear  stress  t  on  the  plane  x=0  corresponding  to  the  sth  mode 
xy 
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/s\  vs 

'  =  G  -g^~  =  -ikgG  cos(Ygz)  exp(iu)t) 


(A3. 7) 


Because  of  symmetry  about  x=0  it  is  sufficient  to  consider  only 
the  region  x>0.  The  shear  stress  T  at  x=0  vanishes  except  at  the 
point  (0,0)  where  the  load  P  =  1  •  exp(itot)  is  applied  as  a  Dirac  delta 
function  S(z)  defined  by 


6(z)  =0  for  z/o 

rZ  y 

/  6(z)dz  =  —  for  all  e>0 

J r\ 


(A3. 8) 


The  complete  solution  is  some  linear  combination  of  the  eigenfunctions 
given  by  Eq.  (A3. 4).  The  boundary  condition  at  x=0  is  therefore 


T  =  'S  a  =  -6(z)  exp(iwt) 

xy  Z-.  s  xy  v 


(A3.?) 


which  by  (A3. 7)  yields 


iG  y  a  k  cos(y  z)  =  6(z) 
Z_/  s  s  s 


(A3. 10) 


The  coefficients  ag  can  be  obtained  by  multiplication  of  Eq.  (A3. 10)  with 
cos(Yrz)  and  integration  over  z 


“  r  „H 

iG  )  a  k  J 

,  L  s  s  o 


cos(ysz)cos(y 


\ 

z)dz  =  / 

JO 


5(z)cos(y^z) dz 


(A3. 11) 


The  integral  on  the  right  has  the  value  1/2  and  the  integrals  in  the  sum 
vanish  except  for  s=r  in  which  case  the  value  is  H/2.  Hence  Eq.  (A3. 11) 


reduces  to 


-i 

as  =  GHk 


s-1 • • • 


(A3 *12) 


which  defines  the  coefficients  ag  and  the  complete  solution  for  x>0  is 


-j  j iwjiiipwpiiBP»iwi!u . ■  ■< 1  ■■*  w’sw^x«\'  vii  i}r^w^>/y^j^j|»^wBWT!iy,*?giB»i8v<By^vy>i^JJiM?wy5y»?<yi9g)^By. 


APPENDIX  4 

COMPUTATION  OF  HANKEL  FUNCTIONS 


It  is  not  advisable  to  compute  values  of  the  Hankel  functions  of 

(2) 

the  second  kind  and  mth  order,  H  (z),  for  complex  arguments, 

m 

z  =  |  z  |  exp(i$)  ,  in  the  range  0<— <*><tt  from  the  values  of  the  Bessel 

functions  J  (z)  and  Y  (z)  by 
m  m 

H(2)(z)  =  J  (z)  -  i  Y  (z)  (A4.1) 

vn  m  m 

because  severe  cancellations  occur  in  the  subtraction  unless  |z|  is 
small  (e.g.  9  significart  digits  would  be  lost  for  |z|  =  10  and 
9  =  -5tt/6)  [2]. 

The  FORTRAN  IV  subroutine  HANKEL,  which  is  listed  in  Appendix  5, 
computes  H^(z)  ar- !  H^(z)  for  0<-4><ir  by  ascending  series  if  | z | <10 
and  by  asymptotic  expansions  if  |z|>10.  The  computations  are  based  on 
the  following  expressions  which  are  taken  from  [2]  and  are  rewritten. 
Ascending  series: 


H^2)  (z)  =  a  +  {a  +  i2/ir  •  (1  +  1/2  +  ...  +l/k)}  bj 
k=l 

H^2) (z)  =  [a  +  i/~  -  (1  +  4/z2)]  z/2 

tJO 

+  z/2  •  ^{a  +  12/tt  •  (1  +  1/2  +  ...+  1/k  +  0.5/(kFl)  ]  }b2/(k+l) 
k=l 

in  which  i  =  Y-l, 

a  =  l-i2/a  •  [0.577215664901533  +  £n(z/2)] 

.  k 

k  -  A5_ 
bk  ~  2k! 


(A4.2) 


(A4.3) 


(A4.4) 

(A4.5) 
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l 
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Asymptotic  expansions: 


H^(z)  =  /2/(ttz)  •  exp(-iz  +  inm/2  +  iir/A)  •  (l  +  ^ 

m  \  ,  -i 

K=1 


(A4.6) 


in  which 


c,  =  Ik!  (i8z)k]_1  *  n  [Am7  -  (2k-l)2] 
k  j=l 


(AA.7) 


In  the  program,  the  series,  Eqs-  (A4.2),  (A4.3),  and  (AA,6),^are 
truncated  when  jb^i  or  JcJ  <10"8-  Computed  values  of  Hq  and 
were  checked  against  tables  133,  3A]  for  several  values  of  (zj  and  <}>. 
The  maximum  relative  error  found  was  about  5  x  10  and  occurred  m  the 
transition  zone,  |zj~10,  for  <5>~-tt/2. 
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APPENDIX  5 
PROGRAM  TORVIB 


Identification: 

Analysis  of  torsional  motions  of  circular  footings  on,  or  embedded 
in,  a  layered,  linearly  elastic  or  viscoelastic  medium. 

Programmed  by  Gunter  Waas,  University  of  California,  Berkeley, 

June  1972. 

Purpose: 

The  program  is  designed  to  analyze  harmonic  torsional  motions  in 
an  axisymmetric  system  which  consists  of  one  or  more  linearly  elastic  or 
viscoelastic  layers  supported  by  a  rough  rigid  base.  The  layered  system 
is  subdivided  into  a  finite  cylindrical  region  I  and  an  infinite  region  R. 
The  two  regions  are  joined  together  at  a  cylindrical  boundary.  The 
harmonic  motions  are  excited  at  nodal  joints  either  by  harmonic  forces 
acting  in  the  angular  direction  0,  i.e.  perpendicular  to  the  axisymmetric 
cross-section,  or  by  prescribed  displacements  in  the  same  direction. 

The  program  consists  of  the  main  program  TORVIB  and  the  9  sub¬ 
routines:  INPUT,  BOUNCO,  STIFFT,  MODIFY,  FORCES,  DISPLA,  ALAMDAB, 

HANKEL  and  CYMSOL. 

Discretization: 

The  region  I  is  discretized  by  finite  elements  which  are  rings 
with  rectangular  cross-sections,  while  the  region  R  is  subdivided  into 
sublayers.  A  typical  mesh  is  presented  in  the  figure  below  which  indicates 
the  order  in  which  elements  and  nodal  joints  are  numbered. 

Size  Limitations: 

The  program  uses  dynamic  storage  allocation  for  array:;,  the  length 
of  which  is  problem  dependent.  The  total  storage  available  for  these 
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arrays  is  defined  by  the  length,  LA,  of  the  array  A  which  appears  in  the 
blank  COMMON  statement  of  program  TORVIB.  The  required  length  for  any 
problem  is 

LA  >  NL2  (2  NC  +  4)  +  NL  (9  NC  +  4  NDI  +  7)  +  9  NDI  +  2  NP  +  NC 


in  which  NL,  NC,  NDI,  and  NP  are  the  numbers  defined  in  paragraph  B 
under  Input  Data. 


Input  Data: 


A.  HEADING  CARD  (12A6)  (Read  in  TORVIB) 

Columns  1  to  72  contain  alphanumeric  data  to  be  printed  as  titles  on 
the  output. 

B.  CONTROL  CARD  (815)  (Read  in  TORVIB) 

Cols-  1-5  No.  of  main  layers,  NML 

6-10  Total  number  of  sublayers,  NL 
11-15  No.  of  columns  of  elements,  NC 

16-20  No.  of  elements,  the  material  properties  of  which  are  to 
be  reset,  NEC 

21-25  No.  of  nodal  loads  (external  forces  or  prescribed  dis¬ 
placements)  ,  NDI 

26-30  No.  of  points  in  Region  R  where  displacements  are  to  be 
computed,  NP 

31-35  Loading  parameter,  MODE 

If  M0DE=1  loading  consists  of  nodal  forces 
If  MQDE=2  loading  consists  of  prescribed  nodal  displace¬ 
ments 

36-40  No.  of  frequencies  at  which  the  system  is  to  be  analyzed, 
NFR 

(E.g.  in  the  figure  NML=3;  NL=5;  NC=3,  NDI=3;  M0DE=2) 

C.  LAYER  CARDS  (I5,4F10.0)  (Read  in  INPUT) 

One  card  for  each  main  layer  commencing  with  the  surface  layer  and 
ending  with  the  bottom  layer. 

Cols.  1-5  No.  of  sublayers  to  the  main  layer,  NS 
6-15  Thickness  of  the  main  layer,  X 
16-25  Mass  density  of  the  main  layer,  Y 
26-35  Real  part  of  the  shear  modulus,  C^ 

36-45  Imaginary  part  of  the  shear  modulus,  C£ 

Note:  These  data  also  define  the  thickness  and  the  material  properties 
of  the  finite  elements  in  Region  I. 

D.  RADIAL  COORDINATE  CARDS  (5F10.0)  (Read  in  INPUT) 

As  many  cards  are  used  as  needed  to  contain  the  radial  coordinates  of 
the  vertical  grid  lines,  commencing  with  the  line  next  to  the  axis 
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of  symmetry  and  ending  with  the  line  which  defines  the  boundary 
between  Region  I  and  Region  R.  (E.g.,  the  radial  coordinates  in 
the  figure  are  those  of  nodal  joints  1,  6,  and  11)  The  number  of 
radial  coordinates,  NC,  is  specified  in  the  Control  Card. 

E.  ELEMENT  MODIFICATION  CARDS  (4I5,3F10.0)  (Read  in  INPUT) 

These  cards  are  omitted  if  NEC=0,  -  Control  Card.  The  material 
properties  of  the  elements  are  air  ./  defined  by  the  Layer  Cards. 
However,  the  material  properties  of  some  elements  may  be  changed  to 
any  other  values  including  zero.  To  do  this  one  card  must  be  input 
for  each  rectangular  block  of  elements  for  which  the  properties  are 
to  be  reset. 

Cols.  1-5  No.  of  uppermost  sublayer  of  block,  IB 
6-10  No.  of  lowermost  sublayer  of  block,  IE 
11-15  No.  of  leftmost  column  of  block,  JB 
16-20  No.  of  rightmost  column  of  block,  JE 
21-30  Mass  density,  RO 
31-40  Real  part  of  shear  modulus, 

41-50  Imaginary  part  of  shear  modulus,  G£ 

(E.g.,  if  the  properties  of  the  elements  13  and  14  are  to  be  changed 
then  IB=3,  IE=4,  JB=3,  JE=3,  and  NEC=2  in  Control  Card. 

F.  LOAD  CARDS  FOR  FORCES  (I5,2F10.0)  (Read  in  INPUT) 

These  cards  are  omitted  if  M0DE=2  in  Control  Card.  One  card  is 
needed  for  each  nodal  joint  at  which  an  external  force  is  applied. 

Cols.  1-5  Number  of  the  nodal  joint  at  which  force  is  applied 
6-15  Real  part  of  force  per  radian 
16-25  Imaginary  part  of  force  per  radian 

G.  LOAD  CARDS  FOR  DISPLACEMENTS  (1015)  (Read  in  INPUT) 

These  cards  are  omitted  if  M0DE=1  in  the  Control  Card.  The  cards 
contain  the  numbers  of  the  nodal  joints  at  which  non-zero  angular 
displacements  are  prescribed.  All  displacements  are  then  auto¬ 
matically  prescribed  to  be  in  phase  and  have  an  amplitude  equal  to 
the  radius  of  the  respective  nodal  joint.  As  many  cards  are  used  as 
needed  to  contain  all  nodal  joints  with  prescribed  non-zero  displace¬ 
ments.  The  number  of  these  nodal  joints  is  equal  to  NDI  specified 
in  the  Control  Card. 

(In  the  example  the  nodal  joints  are  1,  6,  and  11) 

H.  OUTPUT  SPECIFICATION  CARDS  (I5,F10.0)  (Read  in  INPUT) 

Cols.  1-5  Number  of  the  interface  passing  through  the  point  at 
which  the  motion  in  Region  R  is  to  be  computed,  IZ. 
Interfaces  are  numbered  from  the  surface  downward 
commencing  with  1. 

6-15  r-ordinate  of  point  at  which  motion  is  to  be  computed,  RR 
Use  one  card  for  each  point  in  Region  R.  The  number  of  points,  NP, 
is  specified  in  the  Control  Card. 
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|  1.  FREQUENCY  CARDS  (F10.0)  (Read  in  TORVIB) 

f  Cols.  1-10  Angular  frequency  in  radians  per  unit  time 

f  One  card  for  each  frequency  at  which  the  system  is  to  be  analyzed. 

\  The  number  of  frequencies,  NFR,  is  specified  in  the  Control  Card. 

I  Several  problems  may  be  analyzed  in  one  run.  The  data  set  for  each 

F 

\  problem  commences  with  a  Heading  Card. 


Output  Information; 


The  program  prints  the  following  output: 

(i)  Input  and  generated  data. 

(ii)  Nodal  displacements  including  amplitudes  and  phase  angles  for 

all  nodal  joints  in  Region  I  and  for  specified  points  in  Region  R. 

(iii)  Reaction  forces  for  nodal  joints  at  which  non-zero  displace¬ 
ments  are  prescribed,  including  amplitudes  and  phase  angles. 
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APPENDIX  6 
PROGRAM  P LAXLY 


Identification: 

Analysis  of  plane  and  axisynnnetric  layered  media  of  infinite  extent 
subjected  to  harmonic  loads  which  act  within  the  plane  of  the  cross-section. 

Programmed  by  Gunter  Waas,  University  of  California,  Berkeley, 

June  1972. 

Purpose: 

Program  PLAXLY  is  designed  for  the  analysis  of  harmonic  motion  in 
plane  and  axisynnnetric  systems  of  the  type  shown  in  Figs.  4  and  5. 

Materials  may  be  either  linearly  elastic  or  viscoelastic.  The  loads 
(prescribed  nodal  forces  and  nodal  displacements)  vary  harmonically  with 
time.  Static  loads  may  be  analyzed  by  setting  the  frequency  zero.  The 
program  consists  of  a  main  program  PLAXLY  and  15  subroutines:  INPUTD, 
ELSTIF,  QUAD,  BUUMAT,  GENEP ,  SECEVA,  BOMAP,  INVERTC ,  HANKEL,  STIFF, 

BLOCKS,  MODIFY,  OUTPUTD,  RFORCE  and  OUTDIS. 

Discretization: 

The  irregular  region  I,  see  Figs.  4  and  5,  is  subdivided  into  finite 
elements  with  quadrilateral  cross-sections  as  shown  in  Fig.  6.  The 
layered  regions  are  subdivided  into  thin  sublayers  such  that  the  .layer 
interfaces  coincide  with  the  finite  element  nodes  at  the  plane  or 
cylindrical  boundaries  between  the  irregular  and  tie  layered  regions. 

Any  finite  element,  except  those  adjacent  to  the  boundaries  between  the 
irregular  and  layered  regions,  may  have  zero  moduli  and  zero  mass  density, 
while  tne  sublayers  of  the  layered  regions  must  always  have  non-zero 
property  values. 

The  nunJaring  of  the  finite  clement  nodes  is  arbitrary  with  the 
following  restriction.  If  it  exists  the  left  layered  region  L,  see 
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Fig.  4,  is  connected  to  the  first  NUMLL+1  nodes,  where  NUMLL  is  the  number 
of  sublayers  in  the  left  layered  regions.  The  surface  layer  is  connected 

to  nodes  1  and  2,  the  second  layer  to  noder  2  and  3 . and  the  bottom 

layer  (numbered  NUMLL)  to  nodes  NUMLL  and  NUMLL+1.  The  right  layered 
region  R,  if  existing,  is  connected  to  the  last  NUMLR+1  nodes,  where 
NUMLR  is  the  number  of  sublayers  in  the  right  layered  region.  The  sur¬ 
face  layer  is  connected  to  the  nodes  NUMNl’-NUMLL  and  NUMiP-NUMLL+1 ,  ...» 
and  the  bottom  layer  to  NUMJP-1  and  NUMNP,  where  NUMNP  is  the  total 
number  of  nodes.  Because  the  system  is  supported  by  a  rough  rigid  base, 
see  Figs.  4  and  5,  zero  displacements  must  be  specified  for  the  nodes  at 
the  bottom. 

Size  Limitations: 

The  program  uses  dynamic  storage  allocation  for  arrays  the  lengths 
of  which  are  problem  dependent.  The  total  storage  available  for  these 
arrays  is  the  length,  MTOT,  of  the  array  A  which  appears  in  the  blank 
COMMON  state  nent  of  program  PLAXLY.  The  length  requirement  is 

MTOT  >  MBAND  (5  MB  AN  04-7)  +  3  NUMNP  +  2NPBB 

in  which 

MBAND  =  maximum  half-bandwidth  of  global  stiffness  matrix,  i.e. 

twice  the  maximum  difference  between  nodal  point  numbers 
of  nodes  belonging  to  any  element  +2 
NUMNP  =  total  number  of  nodal  points 

NPBB  =  number  of  points  outside  che  finite  element  region  at  which 
motion  is  to  be  computed. 

MTOT  is  set  15000  in  the  version  of  the  program  listed  below.  If  the 
value  of  MTOT  is  reset  then  the  dimension  of  the  array  A  in  the  blank 
COMMON  statement  of  program  PLAXLY  must  also  be  adjusted. 

Output  Information; 

The  program  prints  the  following  output: 

i)  Input  and  generated  data. 

ii)  Nodal  point  displacements  including  amplitudes  and  phase  angles, 

iii)  Reaction  forces  for  nodes  at  which  non-zero  displacements  are 
prescribed,  including  amplitudes  and  phase  angles. 
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A.  START  CARD  (5H)  (Read  in  PLAXLY) 

The  word  START  must  be  punched  in  columns  1  to  5  on  a  separate  card 

at  the  beginning  of  each  problem.  Several  problems  can  be  analyzed 

in  one  run. 

B.  HEADING  CARD  (12A6)  (Read  in  program  PLAXLY) 

Columns  1  to  72  contain  alphanumerical  data  to  be  printed  as  title  on 

the  output. 

C.  CONTROL  CARD  (815)  (Read  in  PLAXLY) 

Cols.  1-5  No.  of  nodal  points,  N’JMNF 

6-10  No.  of  finite  elements,  NUMEL 
11-15  No.  of  different  materials,  NUMMAT 

16-20  No.  of  sublayers  in  Region  L,  NUMLL,  less  or  equal  40 

21-25  No.  cf  sublayers  in  Region  R,  NUMLR,  iess  ox  equal  40 

26-30  No.  of  points  outside  of  Region  I  where  displacements 
are  to  be  computed,  NP3B 

31-35  No.  of  frequencies  at  which  the  model  is  to  be  analyzed, 
NFRQ 

36-40  If  column  40  is  left  blank  problem  is  taken  to  be  axi- 
symmetric  and  NUMLL  must  be  zero.  1  in  column  40 
indicates  a  plane  problem. 

Note:  If  NUMLR  is  input  as  minus  NUMLL  Region  R  is  understood  to  be 

a  mirror  image  of  Region  L,  and  some  computation  can  be  saved. 

D.  MATERIAL  CARDS  (2I5,5F10.0)  (Read  in  INPUTD) 

On.  '-srd  for  each  different  material;  not  more  than  40  materials. 

Cols.  1-5  Material  number,  M 

6-10  Interpretation  parameter,  INTPR 

If  INTPR=0  Cols.  2i-60  contain  moduli 
If  INTPR=1  Cols.  21-60  contain  wave  velocities 
11-20  Mass  density,  RKO 

21-30  Real  component  of  shear  modulus  or  S-wave  velocity,  Gj^ 

31-40  Imaginary  part  cf  shear  modulus  or  S-wave  velocity,  G£ 

41-50  Real  part  of  Poisson's  ratio  or  P-vave  velocity,  NU^ 

51-60  Imaginary  part  of  Poisson’s  ratio  or  P-wave  velocity,  NU2 

E.  IAYER  CARDS  (1015)  (Read  in  INPUTD) 

Cols.  1-5  Material  number  for  sublayer  1 
6-10  Material  number  for  sublayer  2 
etc. 

Note:  This  data  is  grouped  in  sets,  one  set  for  each  layered  region. 

If  there  are  two  layered  regions  the  first  sets  describes 
Region  L  and  the  second  set  Region  R.  If  NUMLL=0  or  NUMLR=0 
the  respective  set  is  omitted.  In  the  axisymmetric  case  NUMLL 
should  be  zero.  Each  set  consists  of  as  many  cards  as  necessary 
to  contain  the  material  numbers  of  the  layers  sequentially, 
starting  with  the  surface  layer  and  ending  with  the  bottom 
layer. 
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F.  NODAL  POINT  CARDS  (2I5.4F10.0.I5)  (Read  in  INPUTD) 

Cols.  1-5  Nodal  point  number,  NL 

6-10  Parameter  indicating  if  displacements  or  forces  are 
specified,  ICODE 

11-20  r-ordinate,  (x-ordinate  in  plane  case),  R(NL) 

21-30  z-ordinate,  Z(NL) 

31-40  U(1,NL)  see  below 
41-50  U(2,NL)  see  below 

51-55  Parameter  for  nodal  point  generation,  INCL 
If  ICODE  is 

0  U(1,NL)  is  specified  force  in  r-direction  and  U(2,NL)  is 
specified  force  in  z-direction 

1  U(1,NL)  is  specified  displacement  in  r-direction  and  U(2,NL)  is 
specified  force  in  z-direction 

2  U(1,NL)  is  specified  force  in  r-direction  and  U(2,NL)  is 
specified  displacement  in  z-direction 

3  U(1,NL)  is  specified  displacement  in  r-direction  and  U(2,NL)  is 
specified  displacement  in  z-direction 

Note:  Nodal  point  cards  need  not  be  input  in  numerical  sequence. 

Suppose  cards  for  nodes  NA  and  NB  are  input  sequentially.  If 
(NA  -  NB)  < 1  then  nodal  point  data  will  be  generated  for  nodes 
NA+INCL,  NA+2INCL,  ....  NB-INCL  where  INCL  is  the  integer 
specified  on  the  card  for  node  NA.  The  coordinates  for  these 
nodal  points  will  be  obtained  by  linear  interpolation  between 
nodes  NA  and  NB.  The  value  of  ICODE  for  the  generated  nodes 
is  0.  If  (NA  -  NB) >  0  no  data  are  generated.  If  left  blank 
INCL  is  taken  as  1. 

G.  OUTPUT  SPECIFICATION  CARDS  (F5.0,F10.0)  (Read  in  INPUTD) 

Cols.  1-5  Number  of  the  interface  passing  through  the  point  at 

which  motion  in  Region  L  or  Region  R  is  to  be  computed, 

RZ.  Interfaces  are  numbered  from  the  surface  downward 
commencing  with  1. 

6-15  Absolute  r-ordinate  (horizontal)  of  point  at  which  motion 
is  to  be  computed,  RZ. 

Use  one  card  for  each  point  in  Region  L  and  Region  R.  Cards  for  points 
in  Region  L  are  placed  first.  The  number  of  points,  NPBB,  is 
specified  by  the  Control  Card,  see  C. 

H.  ELEMENT  CARDS  (715)  (Read  in  EI.STIF) 

Cols.  1-5  Element  number,  MMM  f~~~ - -i^ 

6-10  Nodal  point  i,  IY(1)  /  \ 

11-15  Nodal  point  j,  IY(2) 

16-20  Nodal  point  k,  IY(3) 

21-25  Nodal  point  1,  IY(4)  2 

26-30  Material  identification  number  for  element,  IY(5y.  If 
.left  blank  it  is  taken  as  1. 
oi'-35  Element  generation  parameter,  INCL 
Note:  Ordei  nodal  points  counter  clockwise  when  the  r-axis  points 
horizontally  to  the  right  and  the  z-axis  points  vertically 
downward? . 


rsj*;  ^t^y  ss- 


'kv'i??  jifS^5 


Mote:  Element  data  may  be  generated  automatically  if  the  material 

number  is  the  same  for  all  elements  in  a  series  and  the  nodal 
point  numbers  can  be  obtained  as  follows: 

IYU^  =  IY(1)  .+1NCL 
IY(2)i  =  IY(2)^+INCL 
IY(3)i  =  IY(3)^+INCL 
IYC'.Ji  =  IY  (4)  j+IMCL 

where  j  refers  to  an  element  for  which  the  element  number 
MMMj  is  one  less  than  the  element  number  MMM^.  In  this  case 
only  the  element  card  for  the  first  element  of  the  series  need 
be  input.  However,  the  last  element  (highest  element  number) 
must  always  be  input.  If  INCL  is  left  blank  then  INCL=1. 


I.  FREQUENCY  CARDS  (F10.0)  (Read  in  PLAXLY) 

Cols.  1-10  Frequency  of  the  exciting  forces  and  displacements,  CPS 
Note:  If  CPS  is  positive  it  is  interpreted  as  the  frequency  in  cycles 
per  second.  If  CPS  is  negative  it  is  interpreted  as  the 
angular  frequency  in  radians  per  second.  The  number  of 
frequency  cards  is  equal  to  NFRQ  specified  in  the  Control 
Card,  see  C. 


J.  srop  CARD  (4H)  (Read  in  PLAXLY) 

For  normal  termaintion  of  execution  the  complete  data  deck  (not  each 
individual  problem)  finishes  with  a  card  with  the  word  STOP  punched 
in  columns  1  to  4. 
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